@c Copyright (C) 1996-2022 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

@DOCSTRING(quad)

@DOCSTRING(quad_options)

@DOCSTRING(quadv)

@DOCSTRING(quadl)

@DOCSTRING(quadgk)

@DOCSTRING(quadcc)

@DOCSTRING(integral)

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}.

@DOCSTRING(trapz)

@DOCSTRING(cumtrapz)

@node Orthogonal Collocation
@section Orthogonal Collocation

@DOCSTRING(colloc)

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

@DOCSTRING(dblquad)

@DOCSTRING(triplequad)

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.

@DOCSTRING(quad2d)

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.

@DOCSTRING(integral2)

@DOCSTRING(integral3)

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.