Mercurial > octave-nkf
changeset 12612:16cca721117b stable
doc: Update all documentation for chapter on Numerical Integration
* cumtrapz.m, dblquad.m, quadgk.m, quadl.m, quadv.m, trapz.m,
triplequad.m, quad.cc, quadcc.cc: Improve docstrings.
* Quad-opts.in: Keep code sample together on a single line.
* mk-opts.pl: Update quad-options function description
* octave.texi: Update order of detailmenu to match order in quad.texi.
* quad.txi: Add language about when to use each quad function,
add examples of using trapz.
* aspell-octave.en.pws: Add new spelling words from quad.texi chapter
| author | Rik <octave@nomad.inbox5.com> |
|---|---|
| date | Sun, 17 Apr 2011 19:57:07 -0700 |
| parents | 2846ea58b288 |
| children | 0e79664b969c 3b2e005e4219 |
| files | doc/interpreter/doccheck/aspell-octave.en.pws doc/interpreter/octave.texi doc/interpreter/quad.txi liboctave/Quad-opts.in mk-opts.pl scripts/general/cumtrapz.m scripts/general/dblquad.m scripts/general/quadgk.m scripts/general/quadl.m scripts/general/quadv.m scripts/general/trapz.m scripts/general/triplequad.m src/DLD-FUNCTIONS/quad.cc src/DLD-FUNCTIONS/quadcc.cc |
| diffstat | 14 files changed, 439 insertions(+), 268 deletions(-) [+] |
line wrap: on
line diff
--- a/doc/interpreter/doccheck/aspell-octave.en.pws Thu Apr 14 18:52:54 2011 -0700 +++ b/doc/interpreter/doccheck/aspell-octave.en.pws Sun Apr 17 19:57:07 2011 -0700 @@ -179,6 +179,7 @@ delaunayn DeleteFcn delim +deltaX demi Demmel DeskJet @@ -575,6 +576,7 @@ nocompute nolabel noncommercially +nonsmooth nonzeros noperm normcdf @@ -674,7 +676,11 @@ QQ QRUPDATE qrupdate +quadcc +quadgk +quadl quadpack +quadv Quantile quantile quantiles @@ -806,6 +812,8 @@ Subfunctions subfunctions subinterval +Subintervals +subintervals sublicenses Sublicensing submatrices @@ -866,6 +874,7 @@ tp tpdf traceback +trapz treelayout treeplot tridiagonal
--- a/doc/interpreter/octave.texi Thu Apr 14 18:52:54 2011 -0700 +++ b/doc/interpreter/octave.texi Sun Apr 17 19:57:07 2011 -0700 @@ -657,8 +657,8 @@ Numerical Integration * Functions of One Variable:: +* Orthogonal Collocation:: * Functions of Multiple Variables:: -* Orthogonal Collocation:: Differential Equations
--- a/doc/interpreter/quad.txi Thu Apr 14 18:52:54 2011 -0700 +++ b/doc/interpreter/quad.txi Sun Apr 17 19:57:07 2011 -0700 @@ -20,19 +20,19 @@ @chapter Numerical Integration Octave comes with several built-in functions for computing the integral -of a function numerically. These functions all solve 1-dimensional -integration problems. +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:: -* Orthogonal Collocation:: @end menu @node Functions of One Variable @section Functions of One Variable -Octave supports three different algorithms for computing the integral +Octave supports five different algorithms for computing the integral @tex $$ \int_a^b f(x) d x @@ -45,30 +45,42 @@ @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 Lobatto rule. @item quadgk Numerical integration using an adaptive Gauss-Konrod rule. -@item quadv -Numerical integration using an adaptive vectorized Simpson's rule. - @item quadcc Numerical integration using adaptive Clenshaw-Curtis rules. -@item trapz -Numerical integration using the trapezoidal method. +@item trapz, cumtrapz +Numerical integration of data using the trapezoidal method. @end table @noindent -Besides these functions Octave also allows you to perform cumulative -numerical integration using the trapezoidal method through the -@code{cumtrapz} function. +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. -@DOCSTRING(quad) +@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 (@math{1e^{-6}}--@math{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 -@DOCSTRING(quad_options) Here is an example of using @code{quad} to integrate the function @tex @@ -95,21 +107,22 @@ @example @group function y = f (x) - y = x .* sin (1 ./ x) .* sqrt (abs (1 - 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 call to @code{quad}, but it makes it much easier to generate a -set of points for plotting (because it makes it possible to call the -function with a vector argument to produce a vector result). +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. -Then we simply call quad: +The second step is to call quad with the limits of integration: @example @group -[v, ier, nfun, err] = quad ("f", 0, 3) +[q, ier, nfun, err] = quad ("f", 0, 3) @result{} 1.9819 @result{} 1 @result{} 5061 @@ -121,13 +134,83 @@ 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: + +@example +@group +# Verify integral (x^3) = x^4/4 +f = inline ("x.^3"); +quadgk (f, 0, 1) + @result{} 0.25000 + +# 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(quadv) +@DOCSTRING(quadcc) + +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. -@DOCSTRING(quadcc) +@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) @@ -183,9 +266,9 @@ @section Functions of Multiple Variables Octave does not have built-in functions for computing the integral of -functions of multiple variables directly. It is however possible to +functions of multiple variables directly. It is possible, however, to compute the integral of a function of multiple variables using the -functions for one-dimensional integrals. +existing functions for one-dimensional integrals. To illustrate how the integration can be performed, we will integrate the function @@ -205,23 +288,23 @@ The first approach creates a function that integrates @math{f} with respect to @math{x}, and then integrates that function with respect to -@math{y}. Since @code{quad} is written in Fortran it cannot be called +@math{y}. Because @code{quad} is written in Fortran it cannot be called recursively. This means that @code{quad} cannot integrate a function that calls @code{quad}, and hence cannot be used to perform the double -integration. It is however possible with @code{quadl}, which is what -the following code does. +integration. Any of the other integrators, however, can be used which is +what the following code demonstrates. @example @group -function I = g(y) - I = ones(1, length(y)); - for i = 1:length(y) - f = @@(x) sin(pi.*x.*y(i)).*sqrt(x.*y(i)); - I(i) = quadl(f, 0, 1); +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 = quadl("g", 0, 1) +I = quadgk ("g", 0, 1) @result{} 0.30022 @end group @end example @@ -232,7 +315,7 @@ @example @group -I = dblquad (@@(x, y) sin(pi.*x.*y).*sqrt(x.*y), 0, 1, 0, 1) +I = dblquad (@@(x, y) sin (pi*x.*y) .* sqrt (x.*y), 0, 1, 0, 1) @result{} 0.30022 @end group @end example @@ -241,12 +324,12 @@ @DOCSTRING(triplequad) -The above mentioned approach works but is fairly slow, and that problem -increases exponentially with the dimensionality the problem. Another +The above mentioned approach works, but is fairly slow, and that problem +increases exponentially with the dimensionality of the integral. Another possible solution is to use Orthogonal Collocation as described in the -previous section. 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 +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), @@ -257,13 +340,13 @@ @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=7} points. +following code computes the studied integral using @math{n=8} points. @example @group -f = @@(x,y) sin(pi*x*y').*sqrt(x*y'); -n = 7; -[t, A, B, q] = colloc(n); +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 @@ -272,7 +355,5 @@ @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, a change of variables is needed. +@math{a} and @math{b}, instead of 0 and 1, then a change of variables is needed. - -
--- a/liboctave/Quad-opts.in Thu Apr 14 18:52:54 2011 -0700 +++ b/liboctave/Quad-opts.in Sun Apr 17 19:57:07 2011 -0700 @@ -34,7 +34,7 @@ DOC_ITEM Non-negative relative tolerance. If the absolute tolerance is zero, the relative tolerance must be greater than or equal to -@code{max (50*eps, 0.5e-28)}. +@w{@code{max (50*eps, 0.5e-28)}}. END_DOC_ITEM TYPE = "double" @@ -59,7 +59,7 @@ DOC_ITEM Non-negative relative tolerance for single precision. If the absolute tolerance is zero, the relative tolerance must be greater than or equal to -@code{max (50*eps, 0.5e-28)}. +@w{@code{max (50*eps, 0.5e-28)}}. END_DOC_ITEM TYPE = "float" INIT_VALUE = "::sqrt (FLT_EPSILON)"
--- a/mk-opts.pl Thu Apr 14 18:52:54 2011 -0700 +++ b/mk-opts.pl Sun Apr 17 19:57:07 2011 -0700 @@ -126,6 +126,7 @@ die "mk-opts.pl: unknown command: $_\n" } } + $INCLUDE = "" if not defined $INCLUDE; # Initialize value if required } sub parse_option_block @@ -230,11 +231,12 @@ if (not defined $DOC_STRING) { - $DOC_STRING = "When called with two arguments, this function\\n\\ -allows you set options parameters for the function \@code{$FCN_NAME}.\\n\\ -Given one argument, \@code{$OPT_FCN_NAME} returns the value of the\\n\\ -corresponding option. If no arguments are supplied, the names of all\\n\\ -the available options and their current values are displayed.\\n\\\n"; + $DOC_STRING = "Query or set options for the function \@code{$FCN_NAME}.\\n\\ +When called with no arguments, the names of all available options and\\n\\ +their current values are displayed.\\n\\ +Given one argument, return the value of the corresponding option.\\n\\ +When called with two arguments, \@code{$OPT_FCN_NAME} set the option\\n\\ +\@var{opt} to value \@var{val}."; } } @@ -909,8 +911,11 @@ print <<"_END_EMIT_OPTIONS_FUNCTION_HDR_"; DEFUN_DLD ($OPT_FCN_NAME, args, , "-*- texinfo -*-\\n\\ -\@deftypefn {Loadable Function} {} $OPT_FCN_NAME (\@var{opt}, \@var{val})\\n\\ +\@deftypefn {Loadable Function} {} $OPT_FCN_NAME ()\\n\\ +\@deftypefnx {Loadable Function} {val =} $OPT_FCN_NAME (\@var{opt})\\n\\ +\@deftypefnx {Loadable Function} {} $OPT_FCN_NAME (\@var{opt}, \@var{val})\\n\\ $DOC_STRING\\n\\ +\\n\\ Options include\\n\\ \\n\\ \@table \@code\\n\\
--- a/scripts/general/cumtrapz.m Thu Apr 14 18:52:54 2011 -0700 +++ b/scripts/general/cumtrapz.m Sun Apr 17 19:57:07 2011 -0700 @@ -17,17 +17,25 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {@var{z} =} cumtrapz (@var{y}) -## @deftypefnx {Function File} {@var{z} =} cumtrapz (@var{x}, @var{y}) -## @deftypefnx {Function File} {@var{z} =} cumtrapz (@dots{}, @var{dim}) +## @deftypefn {Function File} {@var{q} =} cumtrapz (@var{y}) +## @deftypefnx {Function File} {@var{q} =} cumtrapz (@var{x}, @var{y}) +## @deftypefnx {Function File} {@var{q} =} cumtrapz (@dots{}, @var{dim}) ## -## Cumulative numerical integration using trapezoidal method. -## @code{cumtrapz (@var{y})} computes the cumulative integral of the -## @var{y} along the first non-singleton dimension. If the argument -## @var{x} is omitted an equally spaced vector is assumed. @code{cumtrapz -## (@var{x}, @var{y})} evaluates the cumulative integral with respect -## to @var{x}. +## 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. ## +## 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). ## @seealso{trapz, cumsum} ## @end deftypefn
--- a/scripts/general/dblquad.m Thu Apr 14 18:52:54 2011 -0700 +++ b/scripts/general/dblquad.m Sun Apr 17 19:57:07 2011 -0700 @@ -17,18 +17,30 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {} dblquad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{tol}, @var{quadf}, @dots{}) -## Numerically evaluate a double integral. The function over with to -## integrate is defined by @code{@var{f}}, and the interval for the -## integration is defined by @code{[@var{xa}, @var{xb}, @var{ya}, -## @var{yb}]}. The function @var{f} must accept a vector @var{x} and a -## scalar @var{y}, and return a vector of the same length as @var{x}. +## @deftypefn {Function File} {} dblquad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}) +## @deftypefnx {Function File} {} dblquad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{tol}) +## @deftypefnx {Function File} {} dblquad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{tol}, @var{quadf}) +## @deftypefnx {Function File} {} 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}. ## -## If defined, @var{tol} defines the absolute tolerance to which to -## which to integrate each sub-integral. +## @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 @math{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{quadgk}. ## ## Additional arguments, are passed directly to @var{f}. To use the default -## value for @var{tol} one may pass an empty matrix. +## value for @var{tol} or @var{quadf} one may pass an empty matrix ([]). ## @seealso{triplequad, quad, quadv, quadl, quadgk, quadcc, trapz} ## @end deftypefn @@ -66,3 +78,4 @@ %!assert (dblquad (@(x,y) exp(-x.^2 - y.^2) , -1, 1, -1, 1, [], @quadgk), pi * erf(1).^2, 1e-6) %!assert (dblquad (@(x,y) exp(-x.^2 - y.^2) , -1, 1, -1, 1, [], @quadl), pi * erf(1).^2, 1e-6) %!assert (dblquad (@(x,y) exp(-x.^2 - y.^2) , -1, 1, -1, 1, [], @quadv), pi * erf(1).^2, 1e-6) +
--- a/scripts/general/quadgk.m Thu Apr 14 18:52:54 2011 -0700 +++ b/scripts/general/quadgk.m Sun Apr 17 19:57:07 2011 -0700 @@ -17,25 +17,29 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {} quadgk (@var{f}, @var{a}, @var{b}, @var{abstol}, @var{trace}) -## @deftypefnx {Function File} {} quadgk (@var{f}, @var{a}, @var{b}, @var{prop}, @var{val}, @dots{}) +## @deftypefn {Function File} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}) +## @deftypefnx {Function File} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}, @var{abstol}) +## @deftypefnx {Function File} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}, @var{abstol}, @var{trace}) +## @deftypefnx {Function File} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}, @var{prop}, @var{val}, @dots{}) ## @deftypefnx {Function File} {[@var{q}, @var{err}] =} quadgk (@dots{}) -## Numerically evaluate integral using adaptive Gauss-Konrod quadrature. +## +## Numerically evaluate the integral of @var{f} from @var{a} to @var{b} +## using adaptive Gauss-Konrod quadrature. +## @var{f} is a function handle, inline function, or string +## containing the name of the function to evaluate. ## The formulation is based on a proposal by L.F. Shampine, ## @cite{"Vectorized adaptive quadrature in @sc{matlab}", Journal of ## Computational and Applied Mathematics, pp131-140, Vol 211, Issue 2, ## Feb 2008} where all function evaluations at an iteration are -## calculated with a single call to @var{f}. Therefore the function -## @var{f} must be of the form @code{@var{f} (@var{x})} and accept -## vector values of @var{x} and return a vector of the same length -## representing the function evaluations at the given values of @var{x}. -## The function @var{f} can be defined in terms of a function handle, -## inline function or string. +## calculated with a single call to @var{f}. Therefore, the function +## @var{f} must be vectorized and must accept a vector of input values @var{x} +## and return an output vector representing the function evaluations at the +## given values of @var{x}. ## -## The bounds of the quadrature @code{[@var{a}, @var{b}]} can be finite -## or infinite and contain weak end singularities. Variable -## transformation will be used to treat infinite intervals and weaken -## the singularities. For example: +## @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) @@ -46,53 +50,57 @@ ## element-by-element operator @code{./} and all user functions to ## @code{quadgk} should do the same. ## -## The absolute tolerance can be passed as a fourth argument in a manner -## compatible with @code{quadv}. Equally the user can request that -## information on the convergence can be printed is the fifth argument -## is logically true. +## The optional argument @var{tol} defines the absolute tolerance used to stop +## the integration procedure. The default value is @math{1e^{-10}}. +## +## 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}. ## -## Alternatively, certain properties of @code{quadgk} can be passed as -## pairs @code{@var{prop}, @var{val}}. Valid properties are +## Alternatively, properties of @code{quadgk} can be passed to the function as +## pairs @code{"@var{prop}", @var{val}}. Valid properties are ## ## @table @code ## @item AbsTol -## Defines the absolute error tolerance for the quadrature. The default +## Define the absolute error tolerance for the quadrature. The default ## absolute tolerance is 1e-10. ## ## @item RelTol -## Defines the relative error tolerance for the quadrature. The default +## Define the relative error tolerance for the quadrature. The default ## relative tolerance is 1e-5. ## ## @item MaxIntervalCount ## @code{quadgk} initially subdivides the interval on which to perform -## the quadrature into 10 intervals. Sub-intervals that have an -## unacceptable error are sub-divided and re-evaluated. If the number of -## sub-intervals exceeds at any point 650 sub-intervals then a poor +## 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 sub-intervals that can exist before exiting. +## number of subintervals that can exist before exiting. ## ## @item WayPoints -## If there exists discontinuities in the first derivative of the -## function to integrate, then these can be flagged with the -## @code{"WayPoints"} property. This forces the ends of a sub-interval -## to fall on the breakpoints of the function and can result in +## Discontinuities in the first derivative of the function to integrate can be +## flagged with the @code{"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, +## computation, or both. For example, ## ## @example -## quadgk (@@(x) abs (1 - x .^ 2), 0, 2, 'Waypoints', 1) +## quadgk (@@(x) abs (1 - x.^2), 0, 2, 'Waypoints', 1) ## @end example ## ## @noindent ## signals the breakpoint in the integrand at @code{@var{x} = 1}. ## ## @item Trace -## If logically true, then @code{quadgk} prints information on the +## If logically true @code{quadgk} prints information on the ## convergence of the quadrature at each iteration. -##@end table +## @end table ## -## If any of @var{a}, @var{b} or @var{waypoints} is complex, then the +## If any of @var{a}, @var{b}, or @var{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, @@ -108,9 +116,10 @@ ## integrates @code{log (z)} along the square defined by @code{[1+1i, ## 1-1i, -1-1i, -1+1i]} ## -## If two output arguments are requested, then @var{err} returns the -## approximate bounds on the error in the integral @code{abs (@var{q} - -## @var{i})}, where @var{i} is the exact value of the integral. +## The result of the integration is returned in @var{q}. +## @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. ## ## @seealso{quad, quadv, quadl, quadcc, trapz, dblquad, triplequad} ## @end deftypefn @@ -281,7 +290,7 @@ 3 .* (b - a) ./ 4 .* (1 - t.^2); endif - ## Split interval into at least 10 sub-interval with a 15 point + ## Split interval into at least 10 subinterval with a 15 point ## Gauss-Kronrod rule giving a minimum of 150 function evaluations while (length (subs) < 11) subs = [subs' ; subs(1:end-1)' + diff(subs') ./ 2, NaN](:)(1 : end - 1); @@ -294,7 +303,7 @@ ## Singularity will cause divide by zero warnings warning ("off", "Octave:divide-by-zero"); - ## Initial evaluation of the integrand on the sub-intervals + ## Initial evaluation of the integrand on the subintervals [q_subs, q_errs] = __quadgk_eval__ (f, subs); q0 = sum (q_subs); err0 = sum (q_errs); @@ -307,8 +316,8 @@ first = true; while (true) - ## Check for sub-intervals that are too small. Test must be - ## performed in untransformed sub-intervals. What is a good + ## Check for subintervals that are too small. Test must be + ## performed in untransformed subintervals. What is a good ## value for this test. Shampine suggests 100*eps if (any (abs (diff (trans (subs), [], 2) / h0) < 100 * myeps)) q = q0; @@ -333,7 +342,7 @@ break; endif - ## Accept the sub-intervals that meet the convergence criteria + ## Accept the subintervals that meet the convergence criteria idx = find (abs (q_errs) < tol .* abs(diff (subs, [], 2)) ./ h); if (first) q = sum (q_subs (idx)); @@ -347,7 +356,7 @@ endif subs(idx,:) = []; - ## If no remaining sub-intervals exit + ## If no remaining subintervals exit if (rows (subs) == 0) break; endif @@ -356,12 +365,12 @@ disp([rows(subs), err, q0]); endif - ## Split remaining sub-intervals in two + ## Split remaining subintervals in two mid = (subs(:,2) + subs(:,1)) ./ 2; subs = [subs(:,1), mid; mid, subs(:,2)]; - ## If the maximum sub-interval count is met accept remaining - ## sub-interval and exit + ## If the maximum subinterval count is met accept remaining + ## subinterval and exit if (rows (subs) > maxint) warning ("quadgk: maximum interval count (%d) met", maxint); q += sum (q_subs); @@ -369,7 +378,7 @@ break; endif - ## Evaluation of the integrand on the remaining sub-intervals + ## Evaluation of the integrand on the remaining subintervals [q_subs, q_errs] = __quadgk_eval__ (f, subs); endwhile
--- a/scripts/general/quadl.m Thu Apr 14 18:52:54 2011 -0700 +++ b/scripts/general/quadl.m Sun Apr 17 19:57:07 2011 -0700 @@ -22,21 +22,28 @@ ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}) ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{}) ## -## Numerically evaluate integral using adaptive Lobatto rule. -## @code{quadl (@var{f}, @var{a}, @var{b})} approximates the integral of -## @code{@var{f}(@var{x})} to machine precision. @var{f} is either a -## function handle, inline function or string containing the name of -## the function to evaluate. The function @var{f} must return a vector -## of output values if given a vector of input values. +## Numerically evaluate the integral of @var{f} from @var{a} to @var{b} +## using an adaptive 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 +## if given a vector of input values. +## +## @var{a} and @var{b} are the lower and upper limits of integration. Both +## limits must be finite. ## -## If defined, @var{tol} defines the relative tolerance to which to -## which to integrate @code{@var{f}(@var{x})}. While if @var{trace} is -## defined, displays the left end point of the current interval, the -## interval length, and the partial integral. +## The optional argument @var{tol} defines the relative tolerance with which +## to perform the integration. The default value is @code{eps}. ## -## Additional arguments @var{p1}, etc., are passed directly to @var{f}. -## To use default values for @var{tol} and @var{trace}, one may pass -## empty matrices. +## The algorithm used by @code{quadl} involves recursively subdividing the +## integration interval. +## If @var{trace} is defined then for each subinterval display: (1) the left +## end of the subinterval, (2) the length of the subinterval, (3) 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 ([]). ## ## Reference: W. Gander and W. Gautschi, @cite{Adaptive Quadrature - ## Revisited}, BIT Vol. 40, No. 1, March 2000, pp. 84--101. @@ -55,7 +62,7 @@ ## * replace global variable terminate2 with local function need_warning ## * add paper ref to docs -function Q = quadl (f, a, b, tol, trace, varargin) +function q = quadl (f, a, b, tol, trace, varargin) need_warning (1); if (nargin < 4) tol = []; @@ -128,7 +135,7 @@ if (is == 0) is = b-a; endif - Q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin{:}); + q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin{:}); endfunction ## ADAPTLOBSTP Recursive function used by QUADL. @@ -141,7 +148,7 @@ ## ## Walter Gautschi, 08/03/98 -function Q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin) +function q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin) h = (b-a)/2; m = (a+b)/2; alpha = sqrt(2/3); @@ -165,12 +172,12 @@ warning ("quadl: required tolerance may not be met"); need_warning (0); endif - Q = i1; + q = i1; if (trace) - disp ([a, b-a, Q]); + disp ([a, b-a, q]); endif else - Q = (adaptlobstp (f, a, mll, fa, fmll, is, trace, varargin{:}) + q = (adaptlobstp (f, a, mll, fa, fmll, is, trace, varargin{:}) + adaptlobstp (f, mll, ml, fmll, fml, is, trace, varargin{:}) + adaptlobstp (f, ml, m, fml, fm, is, trace, varargin{:}) + adaptlobstp (f, m, mr, fm, fmr, is, trace, varargin{:})
--- a/scripts/general/quadv.m Thu Apr 14 18:52:54 2011 -0700 +++ b/scripts/general/quadv.m Sun Apr 17 19:57:07 2011 -0700 @@ -21,32 +21,43 @@ ## @deftypefnx {Function File} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol}) ## @deftypefnx {Function File} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}) ## @deftypefnx {Function File} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{}) -## @deftypefnx {Function File} {[@var{q}, @var{fcnt}] =} quadv (@dots{}) +## @deftypefnx {Function File} {[@var{q}, @var{nfun}] =} quadv (@dots{}) ## ## Numerically evaluate the integral of @var{f} from @var{a} to @var{b} -## using adaptive Simpson's rule. -## @var{f} is either a function handle, inline function or string +## using an adaptive Simpson's rule. +## @var{f} is a function handle, inline function, or string ## containing the name of the function to evaluate. -## The function defined by @var{f} may be a scalar, vector or array-valued. +## @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. ## -## If a value for @var{tol} is given, it defines the tolerance used to stop -## the adaptation procedure, otherwise the default value of 1e-6 is used. +## @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 tolerance used to stop +## the adaptation procedure. The default value is @math{1e^{-6}}. ## -## The algorithm used by @code{quadv}, involves recursively subdividing the -## integration interval and applying Simpson's rule on each sub-interval. -## If @var{trace} is @var{true}, after computing each of these partial -## integrals, display the total number of function evaluations, the left end -## of the sub-interval, the length of the sub-interval and the approximation -## of the integral over the sub-interval. +## 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 @var{f}. -## To use default values for @var{tol} and @var{trace}, one may pass -## empty matrices. +## 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 ([]). ## -##@seealso{quad, quadl, quadgk, quadcc, trapz, dblquad, triplequad} +## The result of the integration is returned in @var{q}. @var{nfun} indicates +## the number of function evaluations that were made. +## +## Note: @code{quadv} is written in Octave's scripting language and can be +## used recursively in @code{dblquad} and @code{triplequad}, unlike the +## similar @code{quad} function. +## @seealso{quad, quadl, quadgk, quadcc, trapz, dblquad, triplequad} ## @end deftypefn -function [q, fcnt] = quadv (f, a, b, tol, trace, varargin) +function [q, nfun] = quadv (f, a, b, tol, trace, varargin) if (nargin < 3) print_usage (); endif @@ -73,7 +84,7 @@ fa = feval (f, a, varargin{:}); fc = feval (f, c, varargin{:}); fb = feval (f, b, varargin{:}); - fcnt = 3; + nfun = 3; ## If have edge singularities, move edge point by eps*(b-a) as ## discussed in Shampine paper used to implement quadgk @@ -87,10 +98,10 @@ h = (b - a); q = (b - a) / 6 * (fa + 4 * fc + fb); - [q, fcnt, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, q, fcnt, abs (h), + [q, nfun, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, q, nfun, abs (h), tol, trace, varargin{:}); - if (fcnt > 10000) + if (nfun > 10000) warning ("maximum iteration count reached"); elseif (isnan (q) || isinf (q)) warning ("infinite or NaN function evaluations were returned"); @@ -99,16 +110,16 @@ endif endfunction -function [q, fcnt, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, q0, - fcnt, hmin, tol, trace, varargin) - if (fcnt > 10000) +function [q, nfun, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, q0, + nfun, hmin, tol, trace, varargin) + if (nfun > 10000) q = q0; else d = (a + c) / 2; e = (c + b) / 2; fd = feval (f, d, varargin{:}); fe = feval (f, e, varargin{:}); - fcnt += 2; + nfun += 2; q1 = (c - a) / 6 * (fa + 4 * fd + fc); q2 = (b - c) / 6 * (fc + 4 * fe + fb); q = q1 + q2; @@ -118,14 +129,14 @@ endif if (trace) - disp ([fcnt, a, b-a, q]); + disp ([nfun, a, b-a, q]); endif ## Force at least one adpative step. - if (fcnt == 5 || abs (q - q0) > tol) - [q1, fcnt, hmin] = simpsonstp (f, a, c, d, fa, fc, fd, q1, fcnt, hmin, + if (nfun == 5 || abs (q - q0) > tol) + [q1, nfun, hmin] = simpsonstp (f, a, c, d, fa, fc, fd, q1, nfun, hmin, tol, trace, varargin{:}); - [q2, fcnt, hmin] = simpsonstp (f, c, b, e, fc, fb, fe, q2, fcnt, hmin, + [q2, nfun, hmin] = simpsonstp (f, c, b, e, fc, fb, fe, q2, nfun, hmin, tol, trace, varargin{:}); q = q1 + q2; endif
--- a/scripts/general/trapz.m Thu Apr 14 18:52:54 2011 -0700 +++ b/scripts/general/trapz.m Sun Apr 17 19:57:07 2011 -0700 @@ -17,15 +17,38 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {@var{z} =} trapz (@var{y}) -## @deftypefnx {Function File} {@var{z} =} trapz (@var{x}, @var{y}) -## @deftypefnx {Function File} {@var{z} =} trapz (@dots{}, @var{dim}) +## @deftypefn {Function File} {@var{q} =} trapz (@var{y}) +## @deftypefnx {Function File} {@var{q} =} trapz (@var{x}, @var{y}) +## @deftypefnx {Function File} {@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. ## -## Numerical integration using trapezoidal method. @code{trapz -## (@var{y})} computes the integral of the @var{y} along the first -## non-singleton dimension. If the argument @var{x} is omitted a -## equally spaced vector is assumed. @code{trapz (@var{x}, @var{y})} -## evaluates the integral with respect to @var{x}. +## 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; +## q = trapz (y) +## @result{} q = 2.525 # No scaling +## q * 0.1 +## @result{} q = 0.2525 # Approximation to integral by scaling +## trapz (x, y) +## @result{} q = 0.2525 # Same result by specifying @var{x} +## @end group +## @end example ## ## @seealso{cumtrapz} ## @end deftypefn
--- a/scripts/general/triplequad.m Thu Apr 14 18:52:54 2011 -0700 +++ b/scripts/general/triplequad.m Sun Apr 17 19:57:07 2011 -0700 @@ -17,23 +17,34 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {} triplequad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb}, @var{tol}, @var{quadf}, @dots{}) -## Numerically evaluate a triple integral. The function over which to -## integrate is defined by @code{@var{f}}, and the interval for the -## integration is defined by @code{[@var{xa}, @var{xb}, @var{ya}, -## @var{yb}, @var{za}, @var{zb}]}. The function @var{f} must accept a -## vector @var{x} and a scalar @var{y}, and return a vector of the same -## length as @var{x}. +## @deftypefn {Function File} {} triplequad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb}) +## @deftypefnx {Function File} {} triplequad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb}, @var{tol}) +## @deftypefnx {Function File} {} triplequad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb}, @var{tol}, @var{quadf}) +## @deftypefnx {Function File} {} 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}. ## -## If defined, @var{tol} defines the absolute tolerance to which to -## which to integrate each sub-integral. +## @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 @math{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{quadgk}. ## ## Additional arguments, are passed directly to @var{f}. To use the default -## value for @var{tol} one may pass an empty matrix. +## value for @var{tol} or @var{quadf} one may pass an empty matrix ([]). ## @seealso{dblquad, quad, quadv, quadl, quadgk, quadcc, trapz} ## @end deftypefn -function Q = triplequad(f, xa, xb, ya, yb, za, zb, tol, quadf, varargin) +function q = triplequad(f, xa, xb, ya, yb, za, zb, tol, quadf, varargin) if (nargin < 7) print_usage (); endif @@ -50,13 +61,13 @@ varargin = {}; endif - Q = dblquad(@(y, z) inner (y, z, f, xa, xb, tol, quadf, varargin{:}),ya, yb, za, zb, tol); + q = dblquad(@(y, z) inner (y, z, f, xa, xb, tol, quadf, varargin{:}),ya, yb, za, zb, tol); endfunction -function Q = __triplequad_inner__ (y, z, f, xa, xb, tol, quadf, varargin) - Q = zeros (size(y)); +function q = __triplequad_inner__ (y, z, f, xa, xb, tol, quadf, varargin) + q = zeros (size(y)); for i = 1 : length (y) - Q(i) = feval (quadf, @(x) f (x, y(i), z, varargin{:}), xa, xb, tol); + q(i) = feval (quadf, @(x) f (x, y(i), z, varargin{:}), xa, xb, tol); endfor endfunction @@ -64,3 +75,4 @@ % !#assert (triplequad (@(x,y,z) exp(-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1, 1, [], @quadgk), pi ^ (3/2) * erf(1).^3, 1e-6) % !#assert (triplequad (@(x,y,z) exp(-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1, 1, [], @quadl), pi ^ (3/2) * erf(1).^3, 1e-6) % !#assert (triplequad (@(x,y,z) exp(-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1, 1, [], @quadv), pi ^ (3/2) * erf(1).^3, 1e-6) +
--- a/src/DLD-FUNCTIONS/quad.cc Thu Apr 14 18:52:54 2011 -0700 +++ b/src/DLD-FUNCTIONS/quad.cc Sun Apr 17 19:57:07 2011 -0700 @@ -174,36 +174,31 @@ DEFUN_DLD (quad, args, nargout, "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{v} =} quad (@var{f}, @var{a}, @var{b})\n\ -@deftypefnx {Loadable Function} {@var{v} =} quad (@var{f}, @var{a}, @var{b}, @var{tol})\n\ -@deftypefnx {Loadable Function} {@var{v} =} quad (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing})\n\ -@deftypefnx {Loadable Function} {[@var{v}, @var{ier}, @var{nfun}, @var{err}] =} quad (@dots{})\n\ -Integrate a nonlinear function of one variable using @sc{quadpack}.\n\ -The first argument is the name of the function, the function handle, or\n\ -the inline function to call to compute the value of the integrand. It\n\ -must have the form\n\ +@deftypefn {Loadable Function} {@var{q} =} quad (@var{f}, @var{a}, @var{b})\n\ +@deftypefnx {Loadable Function} {@var{q} =} quad (@var{f}, @var{a}, @var{b}, @var{tol})\n\ +@deftypefnx {Loadable Function} {@var{q} =} quad (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing})\n\ +@deftypefnx {Loadable Function} {[@var{q}, @var{ier}, @var{nfun}, @var{err}] =} quad (@dots{})\n\ +Numerically evaluate the integral of @var{f} from @var{a} to @var{b} using\n\ +Fortran routines from @w{@sc{quadpack}}. @var{f} is a function handle, inline\n\ +function, or a string containing the name of the function to evaluate.\n\ +The function must have the form @code{y = f (x)} where @var{y} and @var{x}\n\ +are scalars.\n\ \n\ -@example\n\ -y = f (x)\n\ -@end example\n\ -\n\ -@noindent\n\ -where @var{y} and @var{x} are scalars.\n\ -\n\ -The second and third arguments are limits of integration. Either or\n\ -both may be infinite.\n\ +@var{a} and @var{b} are the lower and upper limits of integration. Either\n\ +or both may be infinite.\n\ \n\ The optional argument @var{tol} is a vector that specifies the desired\n\ accuracy of the result. The first element of the vector is the desired\n\ absolute tolerance, and the second element is the desired relative\n\ tolerance. To choose a relative test only, set the absolute\n\ tolerance to zero. To choose an absolute test only, set the relative\n\ -tolerance to zero.\n\ +tolerance to zero. Both tolerances default to @code{sqrt(eps)} or\n\ +approximately @math{1.5e^{-8}}.\n\ \n\ The optional argument @var{sing} is a vector of values at which the\n\ integrand is known to be singular.\n\ \n\ -The result of the integration is returned in @var{v}. @var{ier}\n\ +The result of the integration is returned in @var{q}. @var{ier}\n\ contains an integer error code (0 indicates a successful integration).\n\ @var{nfun} indicates the number of function evaluations that were\n\ made, and @var{err} contains an estimate of the error in the\n\ @@ -212,8 +207,9 @@ The function @code{quad_options} can set other optional\n\ parameters for @code{quad}.\n\ \n\ -Note: because @code{quad} is written in Fortran it\n\ -cannot be called recursively.\n\ +Note: because @code{quad} is written in Fortran it cannot be called\n\ +recursively. This prevents its use in integrating over more than one\n\ +variable by routines @code{dblquad} and @code{triplequad}.\n\ @seealso{quad_options, quadv, quadl, quadgk, quadcc, trapz, dblquad, triplequad}\n\ @end deftypefn") {
--- a/src/DLD-FUNCTIONS/quadcc.cc Thu Apr 14 18:52:54 2011 -0700 +++ b/src/DLD-FUNCTIONS/quadcc.cc Sun Apr 17 19:57:07 2011 -0700 @@ -1479,11 +1479,59 @@ DEFUN_DLD (quadcc, args, nargout, "-*- texinfo -*-\n\ -@deftypefn {Function File} {[@var{int}, @var{err}, @var{nr_points}] =} quadcc (@var{f}, @var{a}, @var{b}, @var{tol})\n\ -@deftypefnx {Function File} {[@var{int}, @var{err}, @var{nr_points}] =} quadcc (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing})\n\ -Numerically evaluate an integral using the doubly-adaptive\n\ -Clenshaw-Curtis quadrature described by P. Gonnet in @cite{Increasing the\n\ -Reliability of Adaptive Quadrature Using Explicit Interpolants}.\n\ +@deftypefn {Function File} {@var{q} =} quadcc (@var{f}, @var{a}, @var{b})\n\ +@deftypefnx {Function File} {@var{q} =} quadcc (@var{f}, @var{a}, @var{b}, @var{tol})\n\ +@deftypefnx {Function File} {@var{q} =} quadcc (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing})\n\ +@deftypefnx {Function File} {[@var{q}, @var{err}, @var{nr_points}] =} quadcc (@dots{})\n\ +Numerically evaluate the integral of @var{f} from @var{a} to @var{b}\n\ +using the doubly-adaptive Clenshaw-Curtis quadrature described by P. Gonnet\n\ +in @cite{Increasing the Reliability of Adaptive Quadrature Using Explicit\n\ +Interpolants}.\n\ +@var{f} is a function handle, inline function, or string\n\ +containing the name of the function to evaluate.\n\ +The function @var{f} must be vectorized and must return a vector of output\n\ +values if given a vector of input values. For example,\n\ +\n\ +@example\n\ + f = @@(x) x .* sin (1./x) .* sqrt (abs (1 - x));\n\ +@end example\n\ +\n\ +@noindent\n\ +which uses the element-by-element `dot' form for all operators.\n\ +\n\ +@var{a} and @var{b} are the lower and upper limits of integration. Either\n\ +or both limits may be infinite. @code{quadcc} handles an inifinite limit\n\ +by substituting the variable of integration with @code{x=tan(pi/2*u)}.\n\ +\n\ +The optional argument @var{tol} defines the relative tolerance used to stop\n\ +the integration procedure. The default value is @math{1e^{-6}}.\n\ +\n\ +The optional argument @var{sing} contains a list of points where the\n\ +integrand has known singularities, or discontinuities\n\ +in any of its derivatives, inside the integration interval.\n\ +For the example above, which has a discontinuity at x=1, the call to\n\ +@code{quadcc} would be as follows\n\ +\n\ +@example\n\ + int = quadcc (f, a, b, 1.0e-6, [ 1 ]);\n\ +@end example\n\ +\n\ +The result of the integration is returned in @var{q}.\n\ +@var{err} is an estimate of the absolute integration error and\n\ +@var{nr_points} is the number of points at which the integrand was evaluated.\n\ +If the adaptive integration did not converge, the value of\n\ +@var{err} will be larger than the requested tolerance. Therefore, it is\n\ +recommended to verify this value for difficult integrands.\n\ +\n\ +@code{quadcc} is capable of dealing with non-numeric\n\ +values of the integrand such as @code{NaN} or @code{Inf}.\n\ +If the integral diverges, and @code{quadcc} detects this, \n\ +then a warning is issued and @code{Inf} or @code{-Inf} is returned.\n\ +\n\ +Note: @code{quadcc} is a general purpose quadrature algorithm\n\ +and, as such, may be less efficient for a smooth or otherwise\n\ +well-behaved integrand than other methods such as @code{quadgk}.\n\ +\n\ The algorithm uses Clenshaw-Curtis quadrature rules of increasing\n\ degree in each interval and bisects the interval if either the\n\ function does not appear to be smooth or a rule of maximum\n\ @@ -1491,57 +1539,6 @@ L2-norm of the difference between two successive interpolations\n\ of the integrand over the nodes of the respective quadrature rules.\n\ \n\ -For example,\n\ -\n\ -@example\n\ - int = quadcc (f, a, b, 1.0e-6);\n\ -@end example\n\ -\n\ -@noindent\n\ -computes the integral of a function @var{f} in the interval\n\ -[@var{a}, @var{b}] to the relative precision of six\n\ -decimal digits.\n\ -The integrand @var{f} should accept a vector argument and return a vector\n\ -result containing the integrand evaluated at each element of the\n\ -argument, for example:\n\ -\n\ -@example\n\ - f = @@(x) x .* sin (1 ./ x) .* sqrt (abs (1 - x));\n\ -@end example\n\ -\n\ -If the integrand has known singularities or discontinuities\n\ -in any of its derivatives inside the interval,\n\ -as does the above example at x=1, these can be specified in\n\ -the additional argument @var{sing} as follows\n\ -\n\ -@example\n\ - int = quadcc (f, a, b, 1.0e-6, [ 1 ]);\n\ -@end example\n\ -\n\ -The two additional output variables @var{err} and @var{nr_points}\n\ -return an estimate of the absolute integration error and\n\ -the number of points at which the integrand was evaluated\n\ -respectively.\n\ -If the adaptive integration did not converge, the value of\n\ -@var{err} will be larger than the requested tolerance. It is\n\ -therefore recommended to verify this value for difficult\n\ -integrands.\n\ -\n\ -If either @var{a} or @var{b} are @code{+/-Inf}, @code{quadcc}\n\ -integrates @var{f} by substituting the variable of integration\n\ -with @code{x=tan(pi/2*u)}.\n\ -\n\ -@code{quadcc} is capable of dealing with non-numeric\n\ -values of the integrand such as @code{NaN} or @code{Inf}\n\ -, as in the above example at x=0.\n\ -If the integral diverges and @code{quadcc} detects this, \n\ -a warning is issued and @code{Inf} or @code{-Inf} is returned.\n\ -\n\ -Note that @code{quadcc} is a general purpose quadrature algorithm\n\ -and as such may be less efficient for smooth or otherwise\n\ -well-behaved integrand than other methods such as\n\ -@code{quadgk} or @code{trapz}.\n\ -\n\ Reference: P. Gonnet, @cite{Increasing the Reliability of Adaptive\n\ Quadrature Using Explicit Interpolants}, ACM Transactions on\n\ Mathematical Software, Vol. 37, Issue 3, Article No. 3, 2010.\n\