@c DO NOT EDIT!  Generated automatically by munge-texi.pl.

@c Copyright (C) 1996-2025 The Octave Project Developers
@c
@c This file is part of Octave.
@c
@c Octave is free software: you can redistribute it and/or modify it
@c under the terms of the GNU General Public License as published by
@c the Free Software Foundation, either version 3 of the License, or
@c (at your option) any later version.
@c
@c Octave is distributed in the hope that it will be useful, but
@c WITHOUT ANY WARRANTY; without even the implied warranty of
@c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
@c GNU General Public License for more details.
@c
@c You should have received a copy of the GNU General Public License
@c along with Octave; see the file COPYING.  If not, see
@c <https://www.gnu.org/licenses/>.

@node Numerical Integration
@chapter 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.

@menu
* Functions of One Variable::
* Orthogonal Collocation::
* Functions of Multiple Variables::
@end menu

@node Functions of One Variable
@section Functions of One Variable

Octave supports five different adaptive quadrature algorithms for computing
the integral
@tex
$$
 \int_a^b f(x) d x
$$
@end tex
of a function @math{f} over the interval from @math{a} to @math{b}.  These are

@table @code
@item quad
Numerical integration based on Gaussian quadrature.

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

@item quadl
Numerical integration using an adaptive @nospell{Lobatto} rule.

@item quadgk
Numerical integration using an adaptive @nospell{Gauss-Konrod} rule.

@item quadcc
Numerical integration using adaptive @nospell{Clenshaw-Curtis} rules.

In addition, the following functions are also provided:

@item integral
A compatibility wrapper function that will choose between @code{quadv} and
@code{quadgk} depending on the integrand and options chosen.

@item trapz, cumtrapz
Numerical integration of data using the trapezoidal method.
@end table

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

@multitable @columnfractions 0.05 0.15 .80
@headitem @tab Function @tab Characteristics
@item @tab quad   @tab Low accuracy with nonsmooth integrands
@item @tab quadv  @tab Medium accuracy with smooth integrands
@item @tab quadl  @tab Medium accuracy with smooth integrands.  Slower than quadgk.
@item @tab quadgk @tab Medium accuracy (1e-6 -- 1e-9) with smooth integrands.
@item @tab        @tab Handles oscillatory functions and infinite bounds
@item @tab quadcc @tab Low to High accuracy with nonsmooth/smooth integrands
@item @tab        @tab Handles oscillatory functions, singularities, and infinite bounds
@end multitable


Here is an example of using @code{quad} to integrate the function
@tex
$$
 f(x) = x \sin (1/x) \sqrt {|1 - x|}
$$
from $x = 0$ to $x = 3$.
@end tex
@ifnottex

@example
  @var{f}(@var{x}) = @var{x} * sin (1/@var{x}) * sqrt (abs (1 - @var{x}))
@end example

@noindent
from @var{x} = 0 to @var{x} = 3.
@end ifnottex

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:

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

Note the use of the `dot' forms of the operators.  This is not necessary for
the @code{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:

@example
@group
[q, ier, nfun, err] = quad ("f", 0, 3)
     @result{} 1.9819
     @result{} 1
     @result{} 5061
     @result{} 1.1522e-07
@end group
@end example

Although @code{quad} returns a nonzero value for @var{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 @qcode{"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:

@example
@group
# Verify gamma function = (n-1)! for n = 4
f = @@(x) x.^3 .* exp (-x);
quadcc (f, 0, Inf)
     @result{} 6.0000
@end group
@end example

@c quad libinterp/corefcn/quad.cc
@anchor{XREFquad}
@html
<span style="display:block; margin-top:-4.5ex;">&nbsp;</span>
@end html


@deftypefn  {} {@var{q} =} quad (@var{f}, @var{a}, @var{b})
@deftypefnx {} {@var{q} =} quad (@var{f}, @var{a}, @var{b}, @var{tol})
@deftypefnx {} {@var{q} =} quad (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing})
@deftypefnx {} {[@var{q}, @var{ier}, @var{nfev}, @var{err}] =} quad (@dots{})
Numerically evaluate the integral of @var{f} from @var{a} to @var{b} using
Fortran routines from @w{@sc{quadpack}}.

@var{f} is a function handle, inline function, or a string containing the
name of the function to evaluate.  The function must have the form @code{y =
f (x)} where @var{y} and @var{x} are scalars.

@var{a} and @var{b} are the lower and upper limits of integration.  Either
or both may be infinite.

The optional argument @var{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 @code{sqrt (eps)} or
approximately 1.5e-8.

The optional argument @var{sing} is a vector of values at which the
integrand is known to be singular.

The result of the integration is returned in @var{q}.

@var{ier} contains an integer error code (0 indicates a successful
integration).

@var{nfev} indicates the number of function evaluations that were
made.

@var{err} contains an estimate of the error in the solution.

The function @code{quad_options} can set other optional parameters for
@code{quad}.

Note: because @code{quad} is written in Fortran it cannot be called
recursively.  This prevents its use in integrating over more than one
variable by routines @code{dblquad} and @code{triplequad}.
@xseealso{@ref{XREFquad_options,,quad_options}, @ref{XREFquadv,,quadv}, @ref{XREFquadl,,quadl}, @ref{XREFquadgk,,quadgk}, @ref{XREFquadcc,,quadcc}, @ref{XREFtrapz,,trapz}, @ref{XREFdblquad,,dblquad}, @ref{XREFtriplequad,,triplequad}}
@end deftypefn


@c quad_options libinterp/corefcn/Quad-opts.cc
@anchor{XREFquad_options}
@html
<span style="display:block; margin-top:-4.5ex;">&nbsp;</span>
@end html


@deftypefn  {} {} quad_options ()
@deftypefnx {} {val =} quad_options (@var{opt})
@deftypefnx {} {} quad_options (@var{opt}, @var{val})
Query or set options for the function @code{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 @var{opt}.

When called with two arguments, @code{quad_options} sets the option
@var{opt} to value @var{val}.

Options include

@table @asis
@item @qcode{"absolute tolerance"}
Absolute tolerance; may be zero for pure relative error test.

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

@item @qcode{"single precision absolute tolerance"}
Absolute tolerance for single precision; may be zero for pure relative
error test.

@item @qcode{"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
@w{@code{max (50*eps, 0.5e-28)}}.
@end table
@end deftypefn


@c quadv scripts/general/quadv.m
@anchor{XREFquadv}
@html
<span style="display:block; margin-top:-4.5ex;">&nbsp;</span>
@end html


@deftypefn  {} {@var{q} =} quadv (@var{f}, @var{a}, @var{b})
@deftypefnx {} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol})
@deftypefnx {} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace})
@deftypefnx {} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{})
@deftypefnx {} {[@var{q}, @var{nfev}] =} quadv (@dots{})

Numerically evaluate the integral of @var{f} from @var{a} to @var{b}
using an adaptive Simpson's rule.

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

@var{a} and @var{b} are the lower and upper limits of integration.  Both
limits must be finite.

The optional argument @var{tol} defines the absolute tolerance used to stop
the adaptation procedure.  The default value is 1e-6.

The algorithm used by @code{quadv} involves recursively subdividing the
integration interval and applying Simpson's rule on each subinterval.
If @var{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 @var{p1}, etc., are passed directly to the function
@var{f}.  To use default values for @var{tol} and @var{trace}, one may pass
empty matrices ([]).

The result of the integration is returned in @var{q}.

The optional output @var{nfev} indicates the total number of function
evaluations performed.

Note: @code{quadv} is written in Octave's scripting language and can be
used recursively in @code{dblquad} and @code{triplequad}, unlike the
@code{quad} function.
@xseealso{@ref{XREFquad,,quad}, @ref{XREFquadl,,quadl}, @ref{XREFquadgk,,quadgk}, @ref{XREFquadcc,,quadcc}, @ref{XREFtrapz,,trapz}, @ref{XREFdblquad,,dblquad}, @ref{XREFtriplequad,,triplequad}, @ref{XREFintegral,,integral}, @ref{XREFintegral2,,integral2}, @ref{XREFintegral3,,integral3}}
@end deftypefn


@c quadl scripts/general/quadl.m
@anchor{XREFquadl}
@html
<span style="display:block; margin-top:-4.5ex;">&nbsp;</span>
@end html


@deftypefn  {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b})
@deftypefnx {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol})
@deftypefnx {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace})
@deftypefnx {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{})
@deftypefnx {} {[@var{q}, @var{nfev}] =} quadl (@dots{})

Numerically evaluate the integral of @var{f} from @var{a} to @var{b} using
an adaptive @nospell{Lobatto} rule.

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

@var{a} and @var{b} are the lower and upper limits of integration.  Both
limits must be finite.

The optional argument @var{tol} defines the absolute tolerance with which
to perform the integration.  The default value is 1e-6.

The algorithm used by @code{quadl} involves recursively subdividing the
integration interval.  If @var{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 @var{p1}, etc., are passed directly to the function
@var{f}.  To use default values for @var{tol} and @var{trace}, one may pass
empty matrices ([]).

The result of the integration is returned in @var{q}.

The optional output @var{nfev} indicates the total number of function
evaluations performed.

Reference: @nospell{W. Gander and W. Gautschi}, @cite{Adaptive Quadrature -
Revisited}, BIT Vol.@: 40, No.@: 1, March 2000, pp.@: 84--101.
@url{https://www.inf.ethz.ch/personal/gander/}
@xseealso{@ref{XREFquad,,quad}, @ref{XREFquadv,,quadv}, @ref{XREFquadgk,,quadgk}, @ref{XREFquadcc,,quadcc}, @ref{XREFtrapz,,trapz}, @ref{XREFdblquad,,dblquad}, @ref{XREFtriplequad,,triplequad}, @ref{XREFintegral,,integral}, @ref{XREFintegral2,,integral2}, @ref{XREFintegral3,,integral3}}
@end deftypefn


@c quadgk scripts/general/quadgk.m
@anchor{XREFquadgk}
@html
<span style="display:block; margin-top:-4.5ex;">&nbsp;</span>
@end html


@deftypefn  {} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b})
@deftypefnx {} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}, @var{abstol})
@deftypefnx {} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}, @var{abstol}, @var{trace})
@deftypefnx {} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}, "@var{prop}", @var{val}, @dots{})
@deftypefnx {} {[@var{q}, @var{err}] =} quadgk (@dots{})

Numerically evaluate the integral of @var{f} from @var{a} to @var{b}
using adaptive @nospell{Gauss-Kronrod} quadrature.

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

@var{a} and @var{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:

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

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

The optional argument @var{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 @code{quadgk} involves subdividing the integration
interval and evaluating each subinterval.  If @var{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 @var{err},
(3) the current estimate for the integral @var{q}.

The behavior of the algorithm can be configured by passing arguments to
@code{quadgk} as pairs @qcode{"@var{prop}", @var{val}}.  Valid properties
are

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

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

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

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

@noindent
will integrate @code{[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 @code{quadgk}.

@item 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
@code{x == 0} and by specifying a waypoint

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

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

If the function has @strong{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
@code{quadgk}.

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

@item MaxIntervalCount
@code{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 @qcode{"MaxIntervalCount"} can be used
to alter the number of subintervals that can exist before exiting.

@item Trace
If logically true @code{quadgk} prints information on the convergence of the
quadrature at each iteration.
@end table

If any of @var{a}, @var{b}, or @var{waypoints} is complex then the
quadrature is treated as a contour integral along a piecewise linear
path defined by
@code{[@var{a}, @var{waypoints}(1), @var{waypoints}(2), @dots{}, @var{b}]}.
In this case the integral is assumed to have no edge singularities.  For
example,

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

@noindent
integrates @code{log (z)} along the square defined by
@code{[1+1i, -1+1i, -1-1i, +1-1i]}.

The result of the integration is returned in @var{q}.

@var{err} is an approximate bound on the error in the integral
@w{@code{abs (@var{q} - @var{I})}}, where @var{I} is the exact value of the
integral.  If the adaptive integration did not converge, the value of
@var{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 @var{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: @nospell{L.F. Shampine},
@cite{"Vectorized adaptive quadrature in @sc{matlab}"}, Journal of
Computational and Applied Mathematics, pp.@: 131--140, Vol 211, Issue 2,
Feb 2008.

@xseealso{@ref{XREFquad,,quad}, @ref{XREFquadv,,quadv}, @ref{XREFquadl,,quadl}, @ref{XREFquadcc,,quadcc}, @ref{XREFtrapz,,trapz}, @ref{XREFdblquad,,dblquad}, @ref{XREFtriplequad,,triplequad}, @ref{XREFintegral,,integral}, @ref{XREFintegral2,,integral2}, @ref{XREFintegral3,,integral3}}
@end deftypefn


@c quadcc libinterp/corefcn/quadcc.cc
@anchor{XREFquadcc}
@html
<span style="display:block; margin-top:-4.5ex;">&nbsp;</span>
@end html


@deftypefn  {} {@var{q} =} quadcc (@var{f}, @var{a}, @var{b})
@deftypefnx {} {@var{q} =} quadcc (@var{f}, @var{a}, @var{b}, @var{tol})
@deftypefnx {} {@var{q} =} quadcc (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing})
@deftypefnx {} {[@var{q}, @var{err}, @var{nr_points}] =} quadcc (@dots{})
Numerically evaluate the integral of @var{f} from @var{a} to @var{b} using
doubly-adaptive @nospell{Clenshaw-Curtis} quadrature.

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

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

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

@var{a} and @var{b} are the lower and upper limits of integration.  Either or
both limits may be infinite.  @code{quadcc} handles an infinite limit by
substituting the variable of integration with @code{x = tan (pi/2*u)}.

The optional argument @var{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 @var{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 @code{quadcc} would be as follows

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

The result of the integration is returned in @var{q}.

@var{err} is an estimate of the absolute integration error.

@var{nr_points} is the number of points at which the integrand was evaluated.

If the adaptive integration did not converge, the value of @var{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 @var{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.

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

Note: @code{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 @code{quadgk}.

The algorithm uses @nospell{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: @nospell{P. Gonnet}, @cite{Increasing the Reliability of Adaptive
Quadrature Using Explicit Interpolants}, @nospell{ACM} Transactions on
Mathematical Software, Vol.@: 37, Issue 3, Article No.@: 3, 2010.
@xseealso{@ref{XREFquad,,quad}, @ref{XREFquadv,,quadv}, @ref{XREFquadl,,quadl}, @ref{XREFquadgk,,quadgk}, @ref{XREFtrapz,,trapz}, @ref{XREFdblquad,,dblquad}, @ref{XREFtriplequad,,triplequad}}
@end deftypefn


@c integral scripts/general/integral.m
@anchor{XREFintegral}
@html
<span style="display:block; margin-top:-4.5ex;">&nbsp;</span>
@end html


@deftypefn  {} {@var{q} =} integral (@var{f}, @var{a}, @var{b})
@deftypefnx {} {@var{q} =} integral (@var{f}, @var{a}, @var{b}, @var{prop}, @var{val}, @dots{})
@deftypefnx {} {[@var{q}, @var{err}] =} integral (@dots{})

Numerically evaluate the integral of @var{f} from @var{a} to @var{b} using
adaptive quadrature.

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

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

@var{a} and @var{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, @code{integral} will perform a straight line
path integral.  Alternatively, a complex domain path can be specified using
the @qcode{"Waypoints"} option (see below).

Additional optional parameters can be specified using
@qcode{"@var{property}", @var{value}} pairs.  Valid properties are:

@table @code
@item Waypoints
Specifies points to be used in defining subintervals of the quadrature
algorithm, or if @var{a}, @var{b}, or @var{waypoints} are complex then
the quadrature is calculated as a contour integral along a piecewise
continuous path.  For more detail, @pxref{XREFquadgk,,@code{quadgk}}.

@item ArrayValued
@code{integral} expects @var{f} to return a scalar value unless
@var{arrayvalued} is specified as true.  This option will cause
@code{integral} to perform the integration over the entire array and return
@var{q} with the same dimensions as returned by @var{f}.  For more detail
@pxref{XREFquadgk,,@code{quadgk}}.

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

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

The optional output @var{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:
@tex
$$error \leq \max \left( AbsTol, RelTol\cdot\vert q\vert \right)$$
@end tex
@ifnottex

@example
@group
  @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|).
@end group
@end example

@end ifnottex

Known @sc{matlab} incompatibilities:

@enumerate
@item
If tolerances are left unspecified, and any integration limits or waypoints
are of type @code{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.
@sc{matlab} leaves the tighter tolerances appropriate for @code{double}
inputs in place regardless of the class of the integration limits.
@end enumerate

@xseealso{@ref{XREFintegral2,,integral2}, @ref{XREFintegral3,,integral3}, @ref{XREFquad,,quad}, @ref{XREFquadgk,,quadgk}, @ref{XREFquadv,,quadv}, @ref{XREFquadl,,quadl}, @ref{XREFquadcc,,quadcc}, @ref{XREFtrapz,,trapz}, @ref{XREFdblquad,,dblquad}, @ref{XREFtriplequad,,triplequad}}
@end deftypefn


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 @code{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).

@example
@group
x = 0:0.1:pi/2;  # Uniformly spaced points
y = cos (x);
trapz (x, y)
     @result{} 0.99666
@end group
@end example

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.

@example
@group
x = 0:0.1:pi/2;  # Uniformly spaced points
y = sin (x);
trapz (x, y)
     @result{} 0.92849
@end group
@end example

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.

@example
@group
x = linspace (0, pi/2, 16);  # Uniformly spaced, but including endpoint
y = sin (x);
trapz (x, y)
     @result{} 0.99909
@end group
@end example

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 @code{trapz}.

@c trapz scripts/general/trapz.m
@anchor{XREFtrapz}
@html
<span style="display:block; margin-top:-4.5ex;">&nbsp;</span>
@end html


@deftypefn  {} {@var{q} =} trapz (@var{y})
@deftypefnx {} {@var{q} =} trapz (@var{x}, @var{y})
@deftypefnx {} {@var{q} =} trapz (@dots{}, @var{dim})

Numerically evaluate the integral of points @var{y} using the trapezoidal
method.

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

If the optional @var{dim} argument is given, operate along this dimension.

Application Note: If @var{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 @math{x^3}
over the range [0, 1] is @math{x^4/4} or 0.25.  The following code uses
@code{trapz} to calculate the integral in three different ways.

@example
@group
x = 0:0.1:1;
y = x.^3;
## No scaling
q = trapz (y)
  @result{} q = 2.5250
## Approximation to integral by scaling
q * 0.1
  @result{} 0.25250
## Same result by specifying @var{x}
trapz (x, y)
  @result{} 0.25250
@end group
@end example

@xseealso{@ref{XREFcumtrapz,,cumtrapz}}
@end deftypefn


@c cumtrapz scripts/general/cumtrapz.m
@anchor{XREFcumtrapz}
@html
<span style="display:block; margin-top:-4.5ex;">&nbsp;</span>
@end html


@deftypefn  {} {@var{q} =} cumtrapz (@var{y})
@deftypefnx {} {@var{q} =} cumtrapz (@var{x}, @var{y})
@deftypefnx {} {@var{q} =} cumtrapz (@dots{}, @var{dim})
Cumulative numerical integration of points @var{y} using the trapezoidal
method.

@w{@code{cumtrapz (@var{y})}}@ computes the cumulative integral of @var{y}
along the first non-singleton dimension.  Where @code{trapz} reports only
the overall integral sum, @code{cumtrapz} reports the current partial sum
value at each point of @var{y}.

When the argument @var{x} is omitted an equally spaced @var{x} vector with
unit spacing (1) is assumed.  @code{cumtrapz (@var{x}, @var{y})} evaluates
the integral with respect to the spacing in @var{x} and the values in
@var{y}.  This is useful if the points in @var{y} have been sampled
unevenly.

If the optional @var{dim} argument is given, operate along this dimension.

Application Note: If @var{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).
@xseealso{@ref{XREFtrapz,,trapz}, @ref{XREFcumsum,,cumsum}}
@end deftypefn


@node Orthogonal Collocation
@section Orthogonal Collocation

@c colloc libinterp/corefcn/colloc.cc
@anchor{XREFcolloc}
@html
<span style="display:block; margin-top:-4.5ex;">&nbsp;</span>
@end html


@deftypefn {} {[@var{r}, @var{amat}, @var{bmat}, @var{q}] =} colloc (@var{n}, "left", "right")
Compute derivative and integral weight matrices for orthogonal collocation.

Reference: @nospell{J. Villadsen}, @nospell{M. L. Michelsen},
@cite{Solution of Differential Equation Models by Polynomial Approximation}.
@end deftypefn


Here is an example of using @code{colloc} to generate weight matrices
for solving the second order differential equation
@tex
$u^\prime - \alpha u^{\prime\prime} = 0$ with the boundary conditions
$u(0) = 0$ and $u(1) = 1$.
@end tex
@ifnottex
@var{u}' - @var{alpha} * @var{u}'' = 0 with the boundary conditions
@var{u}(0) = 0 and @var{u}(1) = 1.
@end ifnottex

First, we can generate the weight matrices for @var{n} points (including
the endpoints of the interval), and incorporate the boundary conditions
in the right hand side (for a specific value of
@tex
$\alpha$).
@end tex
@ifnottex
@var{alpha}).
@end ifnottex

@example
@group
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);
@end group
@end example

Then the solution at the roots @var{r} is

@example
@group
u = [ 0; (at - alpha * bt) \ rhs; 1]
     @result{} [ 0.00; 0.004; 0.01 0.00; 0.12; 0.62; 1.00 ]
@end group
@end example

@node Functions of Multiple Variables
@section 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 @math{f} with respect to @math{x}, and then integrates
that function with respect to @math{y}.  This procedure can be performed
manually using the following example which integrates the function:

@tex
$$
  f(x, y) = \sin(\pi x y)\sqrt{x y}
$$
@end tex
@ifnottex

@example
f(x, y) = sin(pi*x*y) * sqrt(x*y)
@end example

@end ifnottex
for @math{x} and @math{y} between 0 and 1.

Using @code{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 @code{quad} since it is written in Fortran and cannot be
called recursively.)

@example
@group
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)
      @result{} 0.30022
@end group
@end example

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

@example
@group
I = dblquad (@@(x, y) sin (pi*x.*y) .* sqrt (x.*y), 0, 1, 0, 1)
      @result{} 0.30022
@end group
@end example

@c dblquad scripts/general/dblquad.m
@anchor{XREFdblquad}
@html
<span style="display:block; margin-top:-4.5ex;">&nbsp;</span>
@end html


@deftypefn  {} {@var{q} =} dblquad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb})
@deftypefnx {} {@var{q} =} dblquad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{tol})
@deftypefnx {} {@var{q} =} dblquad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{tol}, @var{quadf})
@deftypefnx {} {@var{q} =} dblquad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{tol}, @var{quadf}, @dots{})
Numerically evaluate the double integral of @var{f}.

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

@var{xa}, @var{ya} and @var{xb}, @var{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 @var{tol} defines the absolute tolerance used to
integrate each sub-integral.  The default value is 1e-6.

The optional argument @var{quadf} specifies which underlying integrator
function to use.  Any choice but @code{quad} is available and the default
is @code{quadcc}.

Additional arguments, are passed directly to @var{f}.  To use the default
value for @var{tol} or @var{quadf} one may pass @qcode{':'} or an empty
matrix ([]).
@xseealso{@ref{XREFintegral2,,integral2}, @ref{XREFintegral3,,integral3}, @ref{XREFtriplequad,,triplequad}, @ref{XREFquad,,quad}, @ref{XREFquadv,,quadv}, @ref{XREFquadl,,quadl}, @ref{XREFquadgk,,quadgk}, @ref{XREFquadcc,,quadcc}, @ref{XREFtrapz,,trapz}}
@end deftypefn


@c triplequad scripts/general/triplequad.m
@anchor{XREFtriplequad}
@html
<span style="display:block; margin-top:-4.5ex;">&nbsp;</span>
@end html


@deftypefn  {} {@var{q} =} triplequad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb})
@deftypefnx {} {@var{q} =} triplequad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb}, @var{tol})
@deftypefnx {} {@var{q} =} triplequad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb}, @var{tol}, @var{quadf})
@deftypefnx {} {@var{q} =} triplequad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb}, @var{tol}, @var{quadf}, @dots{})
Numerically evaluate the triple integral of @var{f}.

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

@var{xa}, @var{ya}, @var{za} and @var{xb}, @var{yb}, @var{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 @var{tol} defines the absolute tolerance used to
integrate each sub-integral.  The default value is 1e-6.

The optional argument @var{quadf} specifies which underlying integrator
function to use.  Any choice but @code{quad} is available and the default
is @code{quadcc}.

Additional arguments, are passed directly to @var{f}.  To use the default
value for @var{tol} or @var{quadf} one may pass @qcode{':'} or an empty
matrix ([]).
@xseealso{@ref{XREFintegral3,,integral3}, @ref{XREFintegral2,,integral2}, @ref{XREFdblquad,,dblquad}, @ref{XREFquad,,quad}, @ref{XREFquadv,,quadv}, @ref{XREFquadl,,quadl}, @ref{XREFquadgk,,quadgk}, @ref{XREFquadcc,,quadcc}, @ref{XREFtrapz,,trapz}}
@end deftypefn


The recursive algorithm for quadrature presented above is referred to as
@qcode{"iterated"}.  A separate 2-D integration method is implemented in the
function @code{quad2d}.  This function performs a @qcode{"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.

@c quad2d scripts/general/quad2d.m
@anchor{XREFquad2d}
@html
<span style="display:block; margin-top:-4.5ex;">&nbsp;</span>
@end html


@deftypefn  {} {@var{q} =} quad2d (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb})
@deftypefnx {} {@var{q} =} quad2d (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{prop}, @var{val}, @dots{})
@deftypefnx {} {[@var{q}, @var{err}, @var{iter}] =} quad2d (@dots{})

Numerically evaluate the two-dimensional integral of @var{f} using adaptive
quadrature over the two-dimensional domain defined by @var{xa}, @var{xb},
@var{ya}, @var{yb} using tiled integration.  Additionally, @var{ya} and
@var{yb} may be scalar functions of @var{x}, allowing for the integration
over non-rectangular domains.

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

Additional optional parameters can be specified using
@qcode{"@var{property}", @var{value}} pairs.  Valid properties are:

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

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

@item MaxFunEvals
The maximum number of function calls to the vectorized function @var{f}.
The default value is 5000.

@item Singular
Enable/disable transforms to weaken singularities on the edge of the
integration domain.  The default value is @var{true}.

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

@item FailurePlot
If @code{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 @var{false}.
@end table

Adaptive quadrature is used to minimize the estimate of error until the
following is satisfied:
@tex
$$error \leq \max \left( AbsTol, RelTol\cdot\vert q\vert \right)$$
@end tex
@ifnottex

@example
@group
        @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|)
@end group
@end example

@end ifnottex

The optional output @var{err} is an approximate bound on the error in the
integral @code{abs (@var{q} - @var{I})}, where @var{I} is the exact value
of the integral.  The optional output @var{iter} is the number of vectorized
function calls to the function @var{f} that were used.

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

@example
@group
@var{f} = @@(@var{x},@var{y}) 2*ones (size (@var{x}));
@var{q} = quad2d (@var{f}, 0, 1, 0, 1)
  @result{} @var{q} =  2
@end group
@end example

The result is a volume, which for this constant-value integrand, is just
@code{@var{Length} * @var{Width} * @var{Height}}.

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

@example
@group
@var{f} = @@(@var{x},@var{y}) 2*ones (size (@var{x}));
@var{ymax} = @@(@var{x}) 1 - @var{x};
@var{q} = quad2d (@var{f}, 0, 1, 0, @var{ymax})
  @result{} @var{q} =  1
@end group
@end example

The result is a volume, which for this constant-value integrand
@w{@math{@var{f} = 2}}, is the Triangle Area x Height or
@w{@code{1/2 * @var{Base} * @var{Width} * @var{Height}}}.

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

@example
@group
@var{f} = @@(@var{x},@var{y}) sinc (@var{x}) * sinc (@var{y}));
@var{q} = quad2d (@var{f}, -1, 1, -1, 1)
  @result{} @var{q} =  12.328  (incorrect)
@var{q} = quad2d (@var{f}, -1, 1, -1, 1, "Vectorized", false)
  @result{} @var{q} =  1.390 (correct)
@var{f} = @@(@var{x},@var{y}) sinc (@var{x}) .* sinc (@var{y});
@var{q} = quad2d (@var{f}, -1, 1, -1, 1)
  @result{} @var{q} =  1.390  (correct)
@end group
@end example

The first result is incorrect as the non-elementwise operator between the
sinc functions in @var{f} create unintended matrix multiplications between
the internal integration arrays used by @code{quad2d}.  In the second
result, setting @qcode{"Vectorized"} to false forces @code{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 @var{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 @sc{matlab} incompatibility: If tolerances are left unspecified, and
any integration limits are of type @code{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.  @sc{matlab} leaves the tighter tolerances appropriate
for @code{double} inputs in place regardless of the class of the
integration limits.

Reference: @nospell{L.F. Shampine},
@cite{@sc{matlab} program for quadrature in 2D}, Applied Mathematics and
Computation, pp.@: 266--274, Vol 1, 2008.

@xseealso{@ref{XREFintegral2,,integral2}, @ref{XREFdblquad,,dblquad}, @ref{XREFintegral,,integral}, @ref{XREFquad,,quad}, @ref{XREFquadgk,,quadgk}, @ref{XREFquadv,,quadv}, @ref{XREFquadl,,quadl}, @ref{XREFquadcc,,quadcc}, @ref{XREFtrapz,,trapz}, @ref{XREFintegral3,,integral3}, @ref{XREFtriplequad,,triplequad}}
@end deftypefn


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

@c integral2 scripts/general/integral2.m
@anchor{XREFintegral2}
@html
<span style="display:block; margin-top:-4.5ex;">&nbsp;</span>
@end html


@deftypefn  {} {@var{q} =} integral2 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb})
@deftypefnx {} {@var{q} =} integral2 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{prop}, @var{val}, @dots{})
@deftypefnx {} {[@var{q}, @var{err}] =} integral2 (@dots{})

Numerically evaluate the two-dimensional integral of @var{f} using adaptive
quadrature over the two-dimensional domain defined by @var{xa}, @var{xb},
@var{ya}, @var{yb} (scalars may be finite or infinite).  Additionally,
@var{ya} and @var{yb} may be scalar functions of @var{x}, allowing for
integration over non-rectangular domains.

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

Additional optional parameters can be specified using
@qcode{"@var{property}", @var{value}} pairs.  Valid properties are:

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

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

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

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

Adaptive quadrature is used to minimize the estimate of error until the
following is satisfied:
@tex
$$error \leq \max \left( AbsTol, RelTol\cdot\vert q\vert \right)$$
@end tex
@ifnottex

@example
@group
        @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|)
@end group
@end example

@end ifnottex

@var{err} is an approximate bound on the error in the integral
@code{abs (@var{q} - @var{I})}, where @var{I} is the exact value of the
integral.

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

@example
@group
@var{f} = @@(@var{x},@var{y}) 2*ones (size (@var{x}));
@var{q} = integral2 (@var{f}, 0, 1, 0, 1)
  @result{} @var{q} =  2
@end group
@end example

The result is a volume, which for this constant-value integrand, is just
@w{@code{@var{Length} * @var{Width} * @var{Height}}}.

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

@example
@group
@var{f} = @@(@var{x},@var{y}) 2*ones (size (@var{x}));
@var{ymax} = @@(@var{x}) 1 - @var{x};
@var{q} = integral2 (@var{f}, 0, 1, 0, @var{ymax})
  @result{} @var{q} =  1
@end group
@end example

The result is a volume, which for this constant-value integrand
@w{@math{@var{f} = 2}}, is the Triangle Area x Height or
@w{@code{1/2 * @var{Base} * @var{Width} * @var{Height}}}.

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

@example
@group
@var{f} = @@(@var{x},@var{y}) sinc (@var{x}) * sinc (@var{y}));
@var{q} = integral2 (@var{f}, -1, 1, -1, 1)
  @result{} @var{q} =  12.328  (incorrect)
@var{q} = integral2 (@var{f}, -1, 1, -1, 1, "Vectorized", false)
  @result{} @var{q} =  1.390 (correct)
@var{f} = @@(@var{x},@var{y}) sinc (@var{x}) .* sinc (@var{y});
@var{q} = integral2 (@var{f}, -1, 1, -1, 1)
  @result{} @var{q} =  1.390  (correct)
@end group
@end example

The first result is incorrect as the non-elementwise operator between the
sinc functions in @var{f} create unintended matrix multiplications between
the internal integration arrays used by @code{integral2}.  In the second
result, setting @qcode{"Vectorized"} to false forces @code{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 @var{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 @sc{matlab} incompatibility: If tolerances are left unspecified, and
any integration limits are of type @code{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.  @sc{matlab} leaves the tighter tolerances appropriate
for @code{double} inputs in place regardless of the class of the
integration limits.

Reference: @nospell{L.F. Shampine},
@cite{@sc{matlab} program for quadrature in 2D}, Applied Mathematics and
Computation, pp.@: 266--274, Vol 1, 2008.

@xseealso{@ref{XREFquad2d,,quad2d}, @ref{XREFdblquad,,dblquad}, @ref{XREFintegral,,integral}, @ref{XREFquad,,quad}, @ref{XREFquadgk,,quadgk}, @ref{XREFquadv,,quadv}, @ref{XREFquadl,,quadl}, @ref{XREFquadcc,,quadcc}, @ref{XREFtrapz,,trapz}, @ref{XREFintegral3,,integral3}, @ref{XREFtriplequad,,triplequad}}
@end deftypefn


@c integral3 scripts/general/integral3.m
@anchor{XREFintegral3}
@html
<span style="display:block; margin-top:-4.5ex;">&nbsp;</span>
@end html


@deftypefn  {} {@var{q} =} integral3 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb})
@deftypefnx {} {@var{q} =} integral3 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb}, @var{prop}, @var{val}, @dots{})

Numerically evaluate the three-dimensional integral of @var{f} using
adaptive quadrature over the three-dimensional domain defined by
@var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb} (scalars may
be finite or infinite).  Additionally, @var{ya} and @var{yb} may be
scalar functions of @var{x} and @var{za}, and @var{zb} maybe be scalar
functions of @var{x} and @var{y}, allowing for integration over
non-rectangular domains.

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

Additional optional parameters can be specified using
@qcode{"@var{property}", @var{value}} pairs.  Valid properties are:

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

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

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

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

Adaptive quadrature is used to minimize the estimate of error until the
following is satisfied:
@tex
$$error \leq \max \left( AbsTol, RelTol\cdot\vert q\vert \right)$$
@end tex
@ifnottex

@example
@group
        @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|)
@end group
@end example

@end ifnottex

@var{err} is an approximate bound on the error in the integral
@code{abs (@var{q} - @var{I})}, where @var{I} is the exact value of the
integral.

Example 1 : integrate over a rectangular volume

@example
@group
@var{f} = @@(@var{x},@var{y},@var{z}) ones (size (@var{x}));
@var{q} = integral3 (@var{f}, 0, 1, 0, 1, 0, 1)
  @result{} @var{q} =  1.00000
@end group
@end example

For this constant-value integrand, the result is a volume which is just
@code{@var{Length} * @var{Width} * @var{Height}}.

Example 2 : integrate over a spherical volume

@example
@group
@var{f} = @@(@var{x},@var{y}) ones (size (@var{x}));
@var{ymax} = @@(@var{x}) sqrt (1 - @var{x}.^2);
@var{zmax} = @@(@var{x},@var{y}) sqrt (1 - @var{x}.^2 - @var{y}.^2);
@var{q} = integral3 (@var{f}, 0, 1, 0, @var{ymax}, 0, @var{zmax})
  @result{} @var{q} =  0.52360
@end group
@end example

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

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

@example
@group
@var{f} = @@(@var{x},@var{y}) sinc (@var{x}) * sinc (@var{y}), * sinc (@var{z});
@var{q} = integral3 (@var{f}, -1, 1, -1, 1, -1, 1)
  @result{} @var{q} =  14.535  (incorrect)
@var{q} = integral3 (@var{f}, -1, 1, -1, 1, -1, 1, "Vectorized", false)
  @result{} @var{q} =  1.6388 (correct)
@var{f} = @@(@var{x},@var{y},@var{z}) sinc (@var{x}) .* sinc (@var{y}), .* sinc (@var{z});
@var{q} = integral3 (@var{f}, -1, 1, -1, 1, -1, 1)
  @result{} @var{q} =  1.6388  (correct)
@end group
@end example

The first result is incorrect as the non-elementwise operator between the
sinc functions in @var{f} create unintended matrix multiplications between
the internal integration arrays used by @code{integral3}.  In the second
result, setting @qcode{"Vectorized"} to false forces @code{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 @var{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 @sc{matlab} incompatibility: If tolerances are left unspecified, and
any integration limits are of type @code{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.  @sc{matlab} leaves the tighter tolerances appropriate
for @code{double} inputs in place regardless of the class of the
integration limits.

Reference: @nospell{L.F. Shampine},
@cite{@sc{matlab} program for quadrature in 2D}, Applied Mathematics and
Computation, pp.@: 266--274, Vol 1, 2008.

@xseealso{@ref{XREFtriplequad,,triplequad}, @ref{XREFintegral,,integral}, @ref{XREFquad,,quad}, @ref{XREFquadgk,,quadgk}, @ref{XREFquadv,,quadv}, @ref{XREFquadl,,quadl}, @ref{XREFquadcc,,quadcc}, @ref{XREFtrapz,,trapz}, @ref{XREFintegral2,,integral2}, @ref{XREFquad2d,,quad2d}, @ref{XREFdblquad,,dblquad}}
@end deftypefn


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 (@pxref{Orthogonal Collocation}).  The integral of a
function @math{f(x,y)} for @math{x} and @math{y} between 0 and 1 can be
approximated using @math{n} points by
@tex
$$
 \int_0^1 \int_0^1 f(x,y) d x d y \approx \sum_{i=1}^n \sum_{j=1}^n q_i q_j f(r_i, r_j),
$$
@end tex
@ifnottex
the sum over @code{i=1:n} and @code{j=1:n} of @code{q(i)*q(j)*f(r(i),r(j))},
@end ifnottex
where @math{q} and @math{r} is as returned by @code{colloc (n)}.  The
generalization to more than two variables is straight forward.  The
following code computes the studied integral using @math{n=8} points.

@example
@group
f = @@(x,y) sin (pi*x*y') .* sqrt (x*y');
n = 8;
[t, ~, ~, q] = colloc (n);
I = q'*f(t,t)*q;
      @result{} 0.30022
@end group
@end example

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