## 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
## .
## -*- 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{})
##
## 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).
##
## @code{integral2} is a wrapper for @code{dblquad} 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 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
## absolute tolerance is 1e-10 (1e-5 for single).
## @end table
##
## Known Matlab incompatibilities:
## @table @code
## @item
## 1. @code{integral2} currently only functions over rectangular domains.
## Implementing @var{ya} and @var{yb} as functions of @var{x} is a planned
## future improvement.
## @item
## 2. A @var{'Method'} property is not yet implemented in Octave due to the lack
## of a 'tiled' integration implementation. All integrals are evaluated using an
## equivalent of the 'iterated' method.
## @item
## 3. The underlying 2d integrator only accepts an Absolute Tolerance. As such
## it is not possible to specify @var{'RelTol'}. A default Relative Tolerance
## value of 1e-6 (1e-4 for single) is used unless Relative Tolerance testing is
## disabled by specifying @var{'RelTol'} as @var{'off'}.
##
## @end table
##
## @seealso{quad, quadgk, quadv, quadl, quadcc, trapz, intergral, dblquad,
## triplequad}
## @end deftypefn
function q = integral2 (f, xa, xb, ya, yb, varargin)
## FIXME: it is possible that a non-rectangular domain could be handled by
## overlaying the integrand with a boolean mask function such than
## the integration occurs over a rectangle but regions outside the
## desired domain contribute zero to the integral. This may be an
## inefficient but acceptable hack to get around the rectangular domain
## limit without having to rewrite the integrating function.
## FIXME: implement 'method' property to let the user select between iterated
## and tiled integration. Tiled integration follows the method of
## matlab's quad2d function, currently unimplemented in Octave. Should
## probably just wait for a quad2d implementation to point the
## integral2 wrapper to instead of trying to recreate it here. The
## following can be added to the help docstring once it is functional:
## @item 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.
## FIXME: implement 'reltol' property once there is a good way to pass this
## value to the underlying integrator. The following can be added to
## the help docstring once it is functional:
## @item RelTol
## Define the relative error tolerance for the quadrature. The default
## relative tolerance is 1e-6 (1e-4 for single).
if (nargin < 5 || (mod (nargin, 2) == 0))
print_usage ();
endif
if (! is_function_handle (f))
print_usage ();
endif
if ((! isscalar (xa)) || (! isscalar (xb)) ...
|| (! isscalar (ya)) || (! isscalar (yb)))
print_usage ();
endif
#check for single or double limits so can set appropriate default tolerance
issingle = (isa (xa, "single") || isa (xb, "single") || ...
isa (ya, "single") || 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
intmethod = [];
integfunc = @quadgk;
## check optional parameters to adjust defaults
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 "reltol"
reltol = varargin{idx++};
if (! ischar (reltol))
warning("integral2: RelTol cannot currently be changed from ", ...
" the default. It can be disabled by specifying ", ...\
"'RelTol', 'off'.");
endif
case "abstol"
abstol = varargin{idx++};
if (! isscalar(abstol))
error("integral2: AbsTol value must be a scalar")
endif
case "method"
intmethod = varargin{idx++};
warning (["integral2: alternate integration methods not yet ", ...
"implemented. Method property ignored."]);
otherwise
error ("integral2: unknown property '%s'", prop);
endswitch
endwhile
if (ischar (reltol))
if (strcmp (tolower (reltol), "off"))
#if don't want reltol defined, set integrator to quadcc
integfunc = @quadcc;
else
error ("integral2: unknown RelTol value '%s'", reltol);
endif
endif
q = dblquad (f, xa, xb, ya, yb, abstol, integfunc);
endfunction
%!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', 'off'), 5000, 1e-10);
## tests from dblquad
%! assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, 'AbsTol', 1e-7),
%! 2*log (2), 1e-9);
%! assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, 'RelTol', 'off'),
%! 2*log (2), 1e-10);
%! 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);
## Test input validation
%!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, 'RelTol', 'foo')