## Copyright (C) 2017 Nicholas Jankowski, David Bateman
##
## 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
## .
## -*- texinfo -*-
## @deftypefn {} {@var{q} =} integral2 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb})
## @deftypefnx {} {@var{q} =} integral2 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{prop}, @var{val}, @dots{})
## @deftypefnx {} {[@var{q}, @var{err}] =} integral2 (@dots{})
##
## Numerically evaluate the two-dimensional integral of @var{f} using adaptive
## quadrature over the two-dimensional domain defined by @var{xa}, @var{xb},
## @var{ya}, @var{yb} (scalars may be finite or infinite). Additionally,
## @var{ya} and @var{yb} may be scalar functions of @var{x}, allowing for the
## integration over non-rectangular domains.
##
## @var{f} is a function handle, inline function, or string containing the name
## of the function to evaluate. The function @var{f} must be of 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}.
##
## 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
## value is 1e-10 (1e-5 for single).
##
## @item RelTol
## Define the relative error tolerance for the quadrature. The default
## value is 1e-6 (1e-4 for single).
##
## @item Vectorized
## Option to disable vectorized integration, forcing octave to use only scalar
## inputs when calling the integrand.
## @end table
##
## Adaptive quadrature is used to minimize the estimate of error until the
## following is satisfied:
## @tex
## $$error \leq \max \left( AbsTol, RelTol\cdot\vert q\vert \right)$$
## @end tex
## @ifnottex
##
## @example
## @group
## @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|).
## @end group
## @end example
##
## @end ifnottex
##
## @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.
##
## Known @sc{matlab} incompatibilities:
##
## @enumerate
## @item
## If tolerances are left unspecified, and any integration limits or waypoints
## are of type @code{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.
## @sc{matlab} leaves the tighter tolerances appropriate for @code{double}
## inputs in place regardless of the class of the integration limits.
##
## @item
## @code{integral2} currently does not have an implemented 2d 'tiled'
## integration method. As such the @qcode{Method} property is ignored.
## @end enumerate
##
## @seealso{integral, integral3, quad, quadgk, quadv, quadl, quadcc, trapz,
## dblquad, triplequad}
## @end deftypefn
function [q, err] = integral2 (f, xa, xb, ya, yb, varargin)
if (nargin < 5 || (mod (nargin, 2) == 0))
print_usage ();
endif
if (! is_function_handle (f))
print_usage ();
endif
if (! (isscalar (xa) && isscalar (xb)))
print_usage ();
endif
## Check for single or double limits to set appropriate default tolerance.
issingle = isa ([xa, xb], "single");
issingle = issingle || (!is_function_handle(ya) && isa(ya, "single"));
issingle = issingle || (!is_function_handle(yb) && isa(yb, "single"));
## Set defaults, update with any specified parameters.
if issingle
abstol = 1e-5;
reltol = 1e-4;
else
abstol = 1e-10;
reltol = 1e-6;
endif
vectorized = true;
singular = false;
# leave hooks in place for tiled method, but default all to iterated for now.
method = "iterated";
idx = 1;
while (idx < nargin - 5)
prop = varargin{idx++};
if (! ischar (prop))
error ("integral2: property PROP must be a string");
endif
switch (tolower (prop))
case "abstol"
abstol = varargin{idx++};
if (! (isnumeric (abstol) && isscalar (abstol) && abstol >= 0))
error ("integral2: AbsTol value must be a numeric scalar >= 0");
endif
case "reltol"
reltol = varargin{idx++};
if (! (isnumeric (reltol) && isscalar (reltol) && reltol >= 0))
error ("integral2: RelTol value must be a numeric scalar >= 0");
endif
case "method"
# FIXME: after having a tiled method, remove warning
method = tolower (varargin{idx++});
if strcmp(method, "tiled")
warning (["integral2: tiled method not implemented yet." ...
"Defaulting to iterated method."]);
endif
case "singular"
# Undocumented option to force weakening of edge singularities, but
# only supported by the 'tiled' method
singular = varargin{idx++};
if !islogical (singular)
error ("integral2: 'singular' must be a logical value");
elseif !(strcmp(method, "iterated") || strcmp(method, "auto"))
error ("integral2: method '%s' unrecognized", method)
endif
case "vectorized"
# option to allow unvectorized functions to be treated
vectorized = varargin{idx++};
if !islogical (vectorized)
error ("integral2: 'vectorized' must be a logical value");
endif
otherwise
error ("integral2: unknown property '%s'", prop);
endswitch
endwhile
if ! vectorized
f = @(x, y) arrayfun(f, x, y);
endif
if strcmp(method, "iterated")
if (! (isscalar (ya) && isscalar (yb)))
## FIXME: should this also check whether a ya(x) or yb(x) returns a
## scalar? As it is, a ya = @(x) [x x] returs true to isscalar despite
## having a vector return value
error (["integral2: Non scalar limits currently unsupported with the" ...
"'iterated' method"]);
endif
q = outer_iterated (f, xa, xb, ya, yb, abstol, reltol);
if nargout == 2
warning("integral2: 'iterated' method can not return estimated error");
err = 0;
endif
else
## FIXME: replace with 'tiled' method block. Add the following back into
## the help text:
## Method
## Specifies the two dimensional integration method to be used, with valid
## options being @var{"auto"}, @var{"tiled"}, or @var{"iterated"}.
## @code{integral} will use @var{"auto"} by default, where it will usually
## choose @var{"tiled"} unless any of the integration limits are infinite.
##
## Also add documentaino for 'singular' parameter that currently has no
## effect.
error("integral2: method tiled not implemented yet");
endif
endfunction
function q = outer_iterated (f, xa, xb, ya, yb, abstol, reltol)
# check upper and lower bounds of y
if !is_function_handle(ya)
if isscalar (ya)
ya = @(x) ya * ones(rows(x), columns(x));
else
error ("integral2: 'ya' must be a constant or a (vectorized) function");
endif
endif
if !is_function_handle(yb)
if isscalar (yb)
yb = @(x) yb * ones(rows(x), columns(x));
else
error ("integral2: 'ya' must be a constant or a (vectorized) function");
endif
endif
inner = @inner_iterated;
q = feval (@quadcc, @(x) inner (x, f, ya, yb, abstol, reltol), xa, xb, ...
[abstol, reltol]);
endfunction
function q = inner_iterated (x, f, ya, yb, abstol, reltol)
q = zeros (size (x));
for i = 1 : length (x)
q(i) = feval (@quadcc, @(y) f(x(i), y), ya(x(i)), yb(x(i)), ...
[abstol, reltol]);
endfor
endfunction
% method="auto" which will default to 'tiled'
%!test
%! f = @(x, y) x .* y;
%! assert (integral2 (f, 0, 1, 0, 1), 0.25, 1e-10);
%!test
%! f = @(x, y) 9 * x.^2 + 15 * y.^2;
%! assert (integral2 (f, 0, 5, -5, 0, "AbsTol", 1e-9), 5000, 1e-9);
%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6), 5000, -1e-6);
%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6, "AbsTol", 1e-9),
%! 5000, 1e-9);
## tests from dblquad
%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-7),
%! 2*log (2), 1e-7);
%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "RelTol", 1e-6),
%! 2*log (2), -1e-6);
%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-8,
%! "RelTol", 1e-6), 2*log (2), -1e-6);
%!assert (integral2 (@(x, y) exp (-x.^2 - y.^2) , -1, 1, -1, 1),
%! pi * erf (1).^2, 1e-10);
%!assert (integral2 (@plus, 1, 2, 3, 4), 5, 1e-10);
%!assert (integral2 (@(x,y) 1 ./ (x + y), 0, 1, 0, @(x) 1 - x), 1, -1e-6);
% method="iterated"
%!test
%! f = @(x, y) x .* y;
%! assert (integral2 (f, 0, 1, 0, 1, "method", "iterated"), 0.25, 1e-10);
%!test
%! f = @(x, y) 9 * x.^2 + 15 * y.^2;
%! assert (integral2 (f, 0, 5, -5, 0, "AbsTol", 1e-9, "method", "iterated"),
%! 5000, 1e-9);
%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6, "method", "iterated"),
%! 5000, -1e-6);
%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6, "method", "iterated",
%! "AbsTol", 1e-9),5000, 1e-9);
## tests from dblquad
%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-7, "method",
%! "iterated"), 2*log (2), 1e-7);
%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "RelTol", 1e-6, "method",
%! "iterated"), 2*log (2), -1e-6);
%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-8,
%! "RelTol", 1e-6, "method", "iterated"), 2*log (2), -1e-6);
%!assert (integral2 (@(x, y) exp (-x.^2 - y.^2) , -1, 1, -1, 1, "method",
%! "iterated"), pi * erf (1).^2, 1e-10);
%!assert (integral2 (@plus, 1, 2, 3, 4, "method", "iterated"), 5, 1e-10);
%!assert (integral2 (@(x,y) 1 ./ (x + y), 0, 1, 0, @(x) 1 - x, "method",
%! "iterated"), 1, -1e-6);
## Test input validation
%!error integral2
%!error integral2 (0, 1 ,2 ,3 ,4)
%!error integral2 (@plus)
%!error integral2 (@plus, 1)
%!error integral2 (@plus, 1, 2)
%!error integral2 (@plus, 1, 2, 3)
%!error integral2 (@plus, 1, 2, 3, [4 5])
%!error integral2 (@plus, 1, 2, 3, "test")
%!error integral2 (@plus, 1, 2, 3, 4, "foo")
%!error integral2 (@plus, 1, 2, 3, 4, "foo", "bar")
%!error integral2 (@plus, 1, 2, 3, 4, 99, "bar")
%!error integral2 (@plus, 1, 2, 3, 4, "AbsTol", "foo")
%!error integral2 (@plus, 1, 2, 3, 4, "AbsTol", [1, 2])
%!error integral2 (@plus, 1, 2, 3, 4, "AbsTol", -1)
%!error integral2 (@plus, 1, 2, 3, 4, "RelTol", "foo")
%!error integral2 (@plus, 1, 2, 3, 4, "RelTol", [1, 2])
%!error integral2 (@plus, 1, 2, 3, 4, "RelTol", -1)