Mercurial > octave
changeset 24055:2eae2ad53eb9
Add missing function integral.m (bug #42037).
* integral.m: New function file, primarily a wrapper to quadgk or quadv.
* /scripts/general/module.mk: Add integral.m to build system.
* quad.txi: Add doc reference to integral after description of quad functions.
* NEWS: Annouce new function integral.
* __unimplemented__.m: Remove integral from the list.
author | Nicholas R. Jankowski <jankowskin@asme.org> |
---|---|
date | Tue, 05 Sep 2017 11:25:21 -0400 |
parents | f7710116bf7d |
children | e18bf7459f79 |
files | NEWS doc/interpreter/quad.txi scripts/general/integral.m scripts/general/module.mk scripts/help/__unimplemented__.m |
diffstat | 5 files changed, 221 insertions(+), 7 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Sun Sep 17 15:53:49 2017 -0400 +++ b/NEWS Tue Sep 05 11:25:21 2017 -0400 @@ -48,6 +48,7 @@ gsvd hgtransform humps + integral openvar repelem rticks
--- a/doc/interpreter/quad.txi Sun Sep 17 15:53:49 2017 -0400 +++ b/doc/interpreter/quad.txi Tue Sep 05 11:25:21 2017 -0400 @@ -32,14 +32,14 @@ @node Functions of One Variable @section Functions of One Variable -Octave supports five different algorithms for computing the integral +Octave supports five different adaptive quadrature algorithms for computing +the integral @tex $$ \int_a^b f(x) d x $$ @end tex -of a function @math{f} over the interval from @math{a} to @math{b}. -These are +of a function @math{f} over the interval from @math{a} to @math{b}. These are @table @code @item quad @@ -57,6 +57,12 @@ @item quadcc Numerical integration using adaptive @nospell{Clenshaw-Curtis} rules. +In addition, the following functions are also provided: + +@item integral +A compatibility wrapper function that will choose between @code(quadv) and +@code(quadgk) depending on the integrand and options chosen. + @item trapz, cumtrapz Numerical integration of data using the trapezoidal method. @end table @@ -165,6 +171,8 @@ @DOCSTRING(quadcc) +@DOCSTRING(integral) + Sometimes one does not have the function, but only the raw (x, y) points from which to perform an integration. This can occur when collecting data in an experiment. The @code{trapz} function can integrate these values as shown in
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/general/integral.m Tue Sep 05 11:25:21 2017 -0400 @@ -0,0 +1,208 @@ +## Copyright (C) 2017 Nicholas Jankowski +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {} {@var{q} =} integral (@var{f}, @var{a}, @var{b}) +## @deftypefnx {} {@var{q} =} integral (@var{f}, @var{a}, @var{b}, @var{prop}, @var{val}, @dots{}) +## +## Numerically evaluate the integral of @var{f} from @var{a} to @var{b} using +## adaptive quadrature. +## +## @code{integral} is a wrapper for @code{quadgk} (for scalar integrands) and +## @code{quadv} (for array-valued integrands) intended to provide Matlab +## compatibility. More control of the numerical integration may be achievable +## by calling the various quadrature functions directly. +## +## @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 when given a vector of input values. +## +## @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. If either +## or both limits are complex, @code{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 +## @qcode{"@var{property}", @var{value}} pairs. Valid properties are: +## +## @table @code +## @item AbsTol +## Define the absolute error tolerance for the quadrature. The default +## absolute tolerance is 1e-10 (1e-5 for single). +## +## @item RelTol +## Define the relative error tolerance for the quadrature. The default +## relative tolerance is 1e-6 (1e-4 for single). +## +## @item Waypoints +## Specifies points to be used in defining subintervals of the quadrature +## algorithm, or if @var{a}, @var{b}, or @var{waypoints} are complex then +## the quadrature is calculated as a contour integral along a piecewise +## continuous path. For more detail see @code{quadgk}. +## +## @item ArrayValued +## @code{integral} expects @var{f} to return a scalar value unless +## @var(arrayvalued) is specified as true. This option will cause +## @code{integral} to perform the integration over the entire array and return +## @var{q} with the same dimensions as returned by @var{f}. +## @end table +## +## Implementation Note: As a consequence of using @code{quadgk} and +## @code{quadv}, certain option combinations are currently unsupported. +## @qcode{"ArrayValued"} cannot be combined with @qcode{"RelTol"} or +## @qcode{"Waypoints"}. This is a known incompatibility with Matlab. +## +## @seealso{quad, quadgk, quadv, quadl, quadcc, trapz, dblquad, triplequad} +## @end deftypefn + +function q = integral (f, a, b, varargin) + + if (nargin < 3 || (mod (nargin, 2) == 0)) + print_usage (); + endif + + if (nargin == 3) + ## Pass the simplest case directly to general integrator. + ## Let quadgk function handle input checks on function and limits. + q = quadgk (f, a, b); + else + ## Parse options to determine how to call integrator + abstol = []; + reltol = []; + waypoints = []; + arrayvalued = false; + + idx = 1; + while (idx < nargin - 3) + prop = varargin{idx++}; + if (! ischar (prop)) + error ("integral: property PROP must be a string"); + endif + + switch (tolower (prop)) + case "reltol" + reltol = varargin{idx++}; + case "abstol" + abstol = varargin{idx++}; + case "waypoints" + waypoints = varargin{idx++}(:); + case "arrayvalued" + arrayvalued = varargin{idx++}; + otherwise + error ("integral: unknown property '%s'", prop); + endswitch + endwhile + + if (arrayvalued) + ## FIXME: replace warning with arrayfun(?) call to quadgk + if (! isempty (waypoints)) + warning(["integral: array-valued quadrature routine currently ", ... + "unable to handle WayPoints. WayPoints are ignored."]); + endif + + ## FIXME: remove warning once we have reltol compatible arrayval'd + ## quad or arrayfun call to quadgk. + if (! isempty (reltol)) + warning(["integral: array-valued quadrature only accepts AbsTol.", ... + " RelTol ignored."]); + endif + + ## quadv accepts empty abstol input. + q = quadv (f, a, b, abstol); + + else + ## quadgk will accept empty values of optional parameters + q = quadgk (f, a, b, "AbsTol", abstol, "RelTol", reltol, ... + "WayPoints", waypoints); + + endif + endif + +endfunction + + +## Matlab compatibility tests +%!test +%! f = @(x) exp (-x.^2) .* log (x).^2; +%! emgamma = 0.57721566490153286; +%! exact = (sqrt (pi)*(8*log (2)^2+8*emgamma*log (2)+pi^2+2*emgamma^2))/16; +%! assert (integral (f, 0, Inf), exact, 1e-6); +%! assert (integral (f, 0, Inf, "RelTol", 1e-12), exact, 1e-12); + +%!test # with parameter +%! f = @(x, c) 1 ./ (x.^3 - 2*x - c); +%! assert (integral (@(x) f(x,5), 0, 2), -0.4605015338467329, 1e-12); + +%!test # with tolerances +%! f = @(x) log(x); +%! assert (integral (@(x) f(x), 0, 1, "AbsTol", 1e-6), -1, 1e-6); + +%!test # waypoints +%! f = @(x) 1./(2.*x-1); +%! assert (integral (f, 0, 0, "Waypoints", [1+1i, 1-1i]), -pi*1i, 1e-12); + +%!test # test array function +%! f = @(x) sin ((1:5)*x); +%! assert (integral (f, 0, 1, "ArrayValued", true), 1./[1:5]-cos(1:5)./[1:5], +%! 1e-10); + +%!test +%! f = @(x) x.^5 .* exp (-x) .* sin (x); +%! assert (integral (f, 0, inf, "RelTol", 1e-8, "AbsTol", 1e-12), -15, 2e-12); + +## tests from quadgk +%!test +%!assert (integral (@sin,-pi,pi), 0, 1e-6); +%!assert (integral (inline ("sin"),-pi,pi), 0, 1e-6); +%!assert (integral ("sin",-pi,pi), 0, 1e-6); +%!assert (integral (@sin,-pi,pi), 0, 1e-6); + +%!assert (integral (@sin,-pi,0), -2, 1e-6); +%!assert (integral (@sin,0,pi), 2, 1e-6); +%!assert (integral (@(x) 1./sqrt (x),0,1), 2, 1e-6); +%!assert (integral (@(x) abs (1 - x.^2),0,2, "Waypoints", 1), 2, 1e-6); +%!assert (integral (@(x) 1./(sqrt (x) .* (x+1)),0,Inf), pi, 1e-6); +%!assert (integral (@(z) log (z),1+i,1+i, "WayPoints", [1-i, -1,-i, -1+i]), ... +%! -pi * 1i, 1e-6); +%!assert (integral (@(x) exp (-x .^ 2),-Inf,Inf), sqrt (pi), 1e-6); +%!assert (integral (@(x) exp (-x .^ 2),-Inf,0), sqrt (pi)/2, 1e-6); + +## tests from quadv +## Test vector-valued functions +%!assert (integral (@(x) [(sin (x)), (sin (2*x))], 0, pi, "ArrayValued", 1), ... +%! [2, 0], 1e-5); + +## Test matrix-valued functions +%!test +%! warning ("off", "Octave:divide-by-zero", "local"); +%! assert (integral (@(x) [x, x, x; x, 1./sqrt(x), x; x, x, x], 0, 1, ... +%! "ArrayValued", 1), ... +%! [0.5, 0.5, 0.5; 0.5, 2, 0.5; 0.5, 0.5, 0.5], 1e-5); + +## Test input validation +%!error integral (@sin) +%!error integral (@sin, 0) +%!error integral (@sin, 0, 1, 1e-6, true, 4) +%!error integral (@sin, 0, 1, "DummyArg") +%!error <property PROP must be a string> integral (@sin, 0, 1, 2, 3) +%!error <unknown property 'foo'> integral (@sin, 0, 1, "foo", 3) +%!error integral (@sin, 0, 1, "AbsTol", ones (2,2)) +%!error integral (@sin, 0, 1, "AbsTol", -1) +%!error integral (@sin, 0, 1, "RelTol", ones (2,2)) +%!error integral (@sin, 0, 1, "RelTol", -1)
--- a/scripts/general/module.mk Sun Sep 17 15:53:49 2017 -0400 +++ b/scripts/general/module.mk Tue Sep 05 11:25:21 2017 -0400 @@ -41,6 +41,7 @@ %reldir%/idivide.m \ %reldir%/inputParser.m \ %reldir%/int2str.m \ + %reldir%/integral.m \ %reldir%/interp1.m \ %reldir%/interp2.m \ %reldir%/interp3.m \
--- a/scripts/help/__unimplemented__.m Sun Sep 17 15:53:49 2017 -0400 +++ b/scripts/help/__unimplemented__.m Tue Sep 05 11:25:21 2017 -0400 @@ -59,10 +59,6 @@ txt = ["griddedInterpolant is not implemented. ", ... "Consider using griddata."]; - case "integral" - txt = ["Octave provides many routines for 1-D numerical integration. ", ... - "Consider quadcc, quad, quadv, quadl, quadgk."]; - case "integral2" txt = ["integral2 is not implemented. Consider using dblquad."];