octave: Functions of One Variable
23.1 Functions of One Variable
==============================
Octave supports five different 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.
‘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’.
DONTPRINTYET See also: 
quad_options XREFquad_options, *notequadv:
DONTPRINTYET DONTPRINTYET See also: quad_options XREFquad_options, quadv
XREFquadv, quadl XREFquadl, quadgk XREFquadgk, *noteDONTPRINTYET DONTPRINTYET DONTPRINTYET See also: quad_options XREFquad_options, quadv
XREFquadv, quadl XREFquadl, quadgk XREFquadgk,
quadcc XREFquadcc, trapz XREFtrapz, *notedblquad:
DONTPRINTYET DONTPRINTYET DONTPRINTYET See also: quad_options XREFquad_options, quadv
XREFquadv, quadl XREFquadl, quadgk XREFquadgk,
quadcc XREFquadcc, trapz XREFtrapz, dblquad
XREFdblquad, triplequad XREFtriplequad.
-- : 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.
DONTPRINTYET See also: quad XREFquad, quadl XREFquadl, *noteDONTPRINTYET DONTPRINTYET See also: quad XREFquad, quadl XREFquadl,
quadgk XREFquadgk, quadcc XREFquadcc, *notetrapz:
DONTPRINTYET DONTPRINTYET DONTPRINTYET See also: quad XREFquad, quadl XREFquadl,
quadgk XREFquadgk, quadcc XREFquadcc, trapz
XREFtrapz, dblquad XREFdblquad, *notetriplequad:
DONTPRINTYET DONTPRINTYET DONTPRINTYET See also: quad XREFquad, quadl XREFquadl,
quadgk XREFquadgk, quadcc XREFquadcc, trapz
XREFtrapz, dblquad XREFdblquad, triplequad
XREFtriplequad.
-- : 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.
<http://www.inf.ethz.ch/personal/gander/>
DONTPRINTYET See also: quad XREFquad, quadv XREFquadv, *noteDONTPRINTYET DONTPRINTYET See also: quad XREFquad, quadv XREFquadv,
quadgk XREFquadgk, quadcc XREFquadcc, *notetrapz:
DONTPRINTYET DONTPRINTYET DONTPRINTYET See also: quad XREFquad, quadv XREFquadv,
quadgk XREFquadgk, quadcc XREFquadcc, trapz
XREFtrapz, dblquad XREFdblquad, *notetriplequad:
DONTPRINTYET DONTPRINTYET DONTPRINTYET See also: quad XREFquad, quadv XREFquadv,
quadgk XREFquadgk, quadcc XREFquadcc, trapz
XREFtrapz, dblquad XREFdblquad, triplequad
XREFtriplequad.
-- : 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.
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.
DONTPRINTYET See also: quad XREFquad, quadv XREFquadv, *noteDONTPRINTYET See also: quad XREFquad, quadv XREFquadv,
quadl XREFquadl, quadcc XREFquadcc, trapz XREFtrapz,
dblquad XREFdblquad, triplequad XREFtriplequad.
-- : 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 inifinite limit
by substituting the variable of integration with ‘x = tan
(pi/2*u)’.
The optional argument TOL defines the relative tolerance used to
stop the integration procedure. The default value is 1e-6.
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.0e-6, [ 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.
DONTPRINTYET See also: quad XREFquad, quadv XREFquadv, *noteDONTPRINTYET See also: quad XREFquad, quadv XREFquadv,
quadl XREFquadl, quadgk XREFquadgk, trapz XREFtrapz,
dblquad XREFdblquad, triplequad XREFtriplequad.
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: cumtrapz XREFcumtrapz.
-- : 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: trapz XREFtrapz, cumsum XREFcumsum.
Info Catalog
octave: Numerical Integration
octave: Orthogonal Collocation