(octave.info)Functions of One Variable
23.1 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, a function
handle, or an inline function. These options make it quite easy to do
integration without having to fully define a function in an m-file. For
example:
# Verify integral (x^3) = x^4/4
f = inline ("x.^3");
quadgk (f, 0, 1)
⇒ 0.25000
# 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, NFUN, 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).
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: Note: quad_options, *note quadv:
XREFquadv, Note: quadl, Note: quadgk, Note:
quadcc, Note: trapz, *note dblquad:
XREFdblquad, Note: 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, NFUN] = 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 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: Note: quad, Note: quadl, Note:
quadgk, Note: quadcc, *note trapz:
XREFtrapz, Note: dblquad, *note triplequad:
XREFtriplequad, Note: integral, *note integral2:
XREFintegral2, Note: 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, NFUN] = 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 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: Note: quad, Note: quadv, Note:
quadgk, Note: quadcc, *note trapz:
XREFtrapz, Note: dblquad, *note triplequad:
XREFtriplequad, Note: integral, *note integral2:
XREFintegral2, Note: 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-Konrod 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.
Alternatively, properties of ‘quadgk’ can be passed to the function
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 the above. 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: Note: quad, Note: quadv, Note:
quadl, Note: quadcc, Note: trapz,
Note: dblquad, Note: triplequad, Note:
integral, Note: integral2, Note:
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. 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.
Implementation Note: For Octave versions ≤ 4.2, ‘quadcc’ accepted a
single tolerance argument which specified the relative tolerance.
For versions 4.4 and 5, ‘quadcc’ will issue a warning when called
with a single tolerance argument indicating that the meaning of
this input has changed from relative tolerance to absolute
tolerance. The warning ID for this message is
"Octave:quadcc:RelTol-conversion". The warning may be disabled
with ‘warning ("off", "Octave:quadcc:RelTol-conversion")’.
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: Note: quad, Note: quadv, Note:
quadl, Note: quadgk, Note: trapz,
Note: dblquad, Note: triplequad.
-- : Q = integral (F, A, B)
-- : Q = integral (F, A, B, PROP, VAL, ...)
Numerically evaluate the integral of F from A to B using adaptive
quadrature.
‘integral’ is a wrapper for ‘quadcc’ (general scalar integrands),
‘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:
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.
2. As a consequence of using ‘quadcc’, ‘quadgk’, and ‘quadv’,
certain option combinations are not supported. Currently,
"ArrayValued" cannot be combined with "RelTol" or "Waypoints".
See also: Note: integral2, *note integral3:
XREFintegral3, Note: quad, Note: quadgk,
Note: quadv, Note: quadl, *note quadcc:
XREFquadcc, Note: trapz, Note: dblquad,
Note: 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;
q = trapz (y)
⇒ q = 2.525 # No scaling
q * 0.1
⇒ q = 0.2525 # Approximation to integral by scaling
trapz (x, y)
⇒ q = 0.2525 # Same result by specifying X
See also: Note: 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).
See also: Note: trapz, Note: cumsum.
automatically generated by info2www version 1.2.2.9