Next: Orthogonal Collocation, Up: Numerical Integration [Contents][Index]
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
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).
nfun 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
.
See also: quad_options, quadv, quadl, quadgk, quadcc, trapz, dblquad, triplequad.
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)
.
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 nfun 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.
See also: quad, quadl, quadgk, quadcc, trapz, dblquad, triplequad, integral, integral2, integral3.
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 nfun 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/
See also: quad, quadv, quadgk, quadcc, trapz, dblquad, triplequad, integral, integral2, integral3.
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.
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 tol 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).
MaxIntervalCount
quadgk
initially subdivides the interval on which to perform the
quadrature into 10 intervals. 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.
WayPoints
Discontinuities in the first derivative of the function to integrate can be
flagged with the "WayPoints"
property. This forces the ends of a
subinterval to fall on the breakpoints of the function and can result in
significantly improved estimation of the error in the integral, faster
computation, or both. For example,
quadgk (@(x) abs (1 - x.^2), 0, 2, "Waypoints", 1)
signals the breakpoint in the integrand at x = 1
.
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 continuous
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.
Reference: L.F. Shampine, "Vectorized adaptive quadrature in MATLAB", Journal of Computational and Applied Mathematics, pp. 131–140, Vol 211, Issue 2, Feb 2008.
See also: quad, quadv, quadl, quadcc, trapz, dblquad, triplequad, integral, integral2, integral3.
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. Therefore, it is recommended to verify this value for difficult integrands.
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.
See also: quad, quadv, quadl, quadgk, trapz, dblquad, triplequad.
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), quadgk
(integrals with specified integration
paths), and quadv
(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 quadv
.
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).
Adaptive quadrature is used to minimize the estimate of error until the following is satisfied:
error <= max (AbsTol, RelTol*|q|).
Known MATLAB incompatibilities:
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.
quadcc
, quadgk
, and quadv
,
certain option combinations are not supported. Currently,
"ArrayValued"
cannot be combined with "RelTol"
or
"Waypoints"
.
See also: 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
.
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
See also: cumtrapz.
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).
Next: Orthogonal Collocation, Up: Numerical Integration [Contents][Index]