## Copyright (C) 2017 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} =} quad2d (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb})
## @deftypefnx {} {@var{q} =} quad2d (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{prop}, @var{val}, @dots{})
## @deftypefnx {} {[@var{q}, @var{err}, @var{iter}] =} quad2d (@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} using tiled integration. 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 MaxFunEvals
## The maximum number of function calls to the vectorized function @var{f}.
## The default value is 5000.
##
## @item FailurePlot
## If @code{quad2d} fails to converge to the desired error tolerance before
## @var{MaxFunEvals} is reached, a plot of the areas that still need refinement
## is created. The default value is @var{false}.
##
## @item Singular
## Use transforms to weaken singularities on the edge of the integration
## doamin. Valid options are @var{true} or @var{false}. The default value is
## @var{true}.
##
## @item Vectorized
## Option to disable vectorized integration, forcing octave to use only scalar
## inputs when calling the integrand. The default value is @var{false}.
##
## @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).
## @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
##
## The optional output @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, and @var{iter} is the number of vectorized function calls
## to the function @var{f} that were used.
##
## Known @sc{matlab} incompatibilities:
##
## @enumerate
## @item
## If tolerances are left unspecified, and any integration limits 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.
## @end enumerate
##
## Reference: @nospell{L.F. Shampine},
## @cite{"@sc{matlab} program for quadrature in 2D"}, Applied Mathematics and
## Computation, pp. 266--274, Vol 1, 2008.
##
## @seealso{integral, quad2d, integral3, quad, quadgk, quadv, quadl, quadcc,
## trapz, dblquad, triplequad}
## @end deftypefn
function [q, err, iter] = quad2d (f, xa, xb, ya, yb, varargin)
if ((nargin < 5) || (mod (nargin, 2) == 0))
print_usage ();
endif
## Convert function given as a string to a function handle
if (ischar (f))
f = @(x) feval (f, x);
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 = true;
idx = 1;
maxiter = 5000;
failureplot = false;
while (idx < nargin - 5)
prop = varargin{idx++};
if (! ischar (prop))
error ("quad2d: property PROP must be a string");
endif
switch (tolower (prop))
case "abstol"
abstol = varargin{idx++};
if (! ((isnumeric (abstol)) && (isscalar (abstol)) && (abstol >= 0)))
error ("quad2d: AbsTol value must be a numeric scalar >= 0");
endif
case "reltol"
reltol = varargin{idx++};
if (! ((isnumeric (reltol)) && (isscalar (reltol)) && (reltol >= 0)))
error ("quad2d: RelTol value must be a numeric scalar >= 0");
endif
case "vectorized"
vectorized = varargin{idx++};
if (! islogical (vectorized))
error ("quad2d: 'vectorized' must be a logical value");
endif
case "singular"
singular = varargin{idx++};
if (! islogical (singular))
error ("quad2d: 'singular' must be a logical value");
endif
case "maxfunevals"
maxiter = varargin{idx++};
if (! ((isnumeric (maxiter)) && (isscalar (maxiter)) &&
(floor (maxiter) == maxiter) && (maxiter >= 1)))
error ("quad2d: MaxFunEvals value must be an integer scalar >= 1");
endif
case "failureplot"
failureplot = varargin{idx++};
if (! islogical (failureplot))
error ("quad2d: 'failureplot' must be a logical value");
endif
otherwise
error ("quad2d: unknown property '%s'", prop);
endswitch
endwhile
if (! vectorized)
f = @(x, y) arrayfun (f, x, y);
endif
# 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 ("quad2d: '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 ("quad2d: 'ya' must be a constant or a (vectorized) function.");
endif
endif
iter = 0;
qaccept = 0;
qerraccept = 0;
if (singular)
# Shampine suggests using the singularity weakening transform
# suggested by Havie
# \int_a^b f(x) dx = \int_0^pi f (g(t)) (dx / dt) dt
# where
# g(t) = ((a - b) * cos(t) + (a + b)) / 2
# dx = - (a - b) * sin(t) / 2 dt
# Now our integral is
# \int_a^b \int_0^1 f(x,y) dydx
# as we already subsitute for "y", so
# gx(tx) = ((a - b) * cos(tx) + (a + b)) / 2
# gy(ty) = (1 - cos(ty)) / 2
# dydx = (b - a) * sin(tx) * sin(ty) / 4 dtydtx
xtrans = @(tx) ((xa - xb) .* cos (tx) + (xa + xb)) ./ 2;
ytrans = @(ty) (1 - cos (ty)) ./ 2;
ztrans = @(tx, ty) (xb - xa) .* sin(tx) .* sin(ty) ./ 4;
area = pi ^ 2;
# Initialize tile list
tilelist(1) = struct("xa", 0, "xb", pi, "ya", 0, "yb", pi, "q", 0, ...
"qerr", Inf);
else
xtrans = @(tx) tx;
ytrans = @(ty) ty;
ztrans = @(tx, ty) 1;
area = (xb - xa);
# Initialize tile list
tilelist(1) = struct("xa", xa, "xb", xb, "ya", 0, "yb", 1, "q", 0, ...
"qerr", Inf);
endif
while ((length (tilelist) > 0) && (iter < maxiter))
# Get tile with the largest error
[~, indx] = max ([tilelist.qerr]);
tile = tilelist(indx);
tilelist(indx) = [];
# Subdivide the tile into 4 subtiles
iter += 4;
tiles(4) = struct("xa", tile.xa, "xb", tile.xa + (tile.xb - tile.xa) / 2,
"ya", tile.ya, "yb", tile.ya + (tile.yb - tile.ya) / 2,
"q", 0, "qerr", 0);
tiles(3) = struct("xa", tile.xa, "xb", tile.xa + (tile.xb - tile.xa) / 2,
"ya", tile.ya + (tile.yb - tile.ya) / 2, "yb", tile.yb,
"q", 0, "qerr", 0);
tiles(2) = struct("xa", tile.xa + (tile.xb - tile.xa) / 2, "xb", tile.xb,
"ya", tile.ya, "yb", tile.ya + (tile.yb - tile.ya) / 2,
"q", 0, "qerr", 0);
tiles(1) = struct("xa", tile.xa + (tile.xb - tile.xa) / 2, "xb", tile.xb,
"ya", tile.ya + (tile.yb - tile.ya) / 2, "yb", tile.yb,
"q", 0, "qerr", 0);
# Perform the quadrature of 4 subtiles
for i = 1:4
[tiles(i).q, tiles(i).qerr] = tensorproduct (f, ya, yb, tiles(i), ...
xtrans, ytrans, ztrans, singular);
endfor
q = qaccept + sum ([[tilelist.q], tiles.q]);
err = qerraccept + sum ([[tilelist.qerr], tiles.qerr]);
tol = max (abstol, reltol .* abs (q));
# Shampine suggests taking a margin of a factor of 8 for
# the local tolerance. That and the fact that we are subdividing
# into 4 tiles, means we divide by 32 at this point.
localtol = tol * ([tile.xb] - [tile.xa]) * ([tile.yb] - [tile.ya]) ...
/ area / 32;
# If global tolerance is met, return
if (err < tol)
break;
endif
# Accept the tiles meeting the tolerance, and add the others back to
# the list of tiles to treat
indx = find ([tiles.qerr] < localtol);
qaccept += sum ([tiles(indx).q]);
qerraccept += sum ([tiles(indx).qerr]);
tiles(indx) = [];
tilelist = [tilelist, tiles];
endwhile
if (iter >= maxiter)
if (err > max (abstol, reltol .* abs (q)))
warning ("quad2d : Maximum number of sub-tiles without convergence.");
else
warning(["quad2d : Maximum number of sub-tiles reached, accuracy ", ...
"may be low."]);
endif
if (failureplot)
newplot;
title("QUAD2D : Areas needing refinement");
for tile = tilelist
xaa = xtrans(tile.xa);
xbb = xtrans(tile.xb);
y1 = ya(xaa) + ytrans(tile.ya) * (yb(xaa) - ya(xaa));
y2 = ya(xaa) + ytrans(tile.yb) * (yb(xaa) - ya(xaa));
y3 = ya(xbb) + ytrans(tile.yb) * (yb(xbb) - ya(xbb));
y4 = ya(xbb) + ytrans(tile.ya) * (yb(xbb) - ya(xbb));
patch ([xaa, xaa, xbb, xbb, xaa],[y1, y2, y3, y4, y1], 'b');
endfor
endif
endif
endfunction
function [q, qerr] = tensorproduct (f, ya, yb, tile, xtrans, ytrans, ...
ztrans, singular)
# The Shampine TwoD paper proposes using a G3,K7 rule in a tensor product
# I couldn't find a tabulated abscissas and weights of a G3,K7 rule publically
# available, so I'm using a G7,K15 rule from Octave's implementation of quadgk
persistent abscissa = [-0.9914553711208126e+00, -0.9491079123427585e+00, ...
-0.8648644233597691e+00, -0.7415311855993944e+00, ...
-0.5860872354676911e+00, -0.4058451513773972e+00, ...
-0.2077849550078985e+00, 0.0000000000000000e+00, ...
0.2077849550078985e+00, 0.4058451513773972e+00, ...
0.5860872354676911e+00, 0.7415311855993944e+00, ...
0.8648644233597691e+00, 0.9491079123427585e+00, ...
0.9914553711208126e+00];
persistent weights15 = ...
[0.2293532201052922e-01, 0.6309209262997855e-01, ...
0.1047900103222502e+00, 0.1406532597155259e+00, ...
0.1690047266392679e+00, 0.1903505780647854e+00, ...
0.2044329400752989e+00, 0.2094821410847278e+00, ...
0.2044329400752989e+00, 0.1903505780647854e+00, ...
0.1690047266392679e+00, 0.1406532597155259e+00, ...
0.1047900103222502e+00, 0.6309209262997855e-01, ...
0.2293532201052922e-01];
persistent weights7 = [0.0, ...
0.1294849661688697e+00, 0.0, ...
0.2797053914892767e+00, 0.0, ...
0.3818300505051889e+00, 0.0, ...
0.4179591836734694e+00, 0.0, ...
0.3818300505051889e+00, 0.0, ...
0.2797053914892767e+00, 0.0, ...
0.1294849661688697e+00, 0.0];
xaa = tile.xa;
xbb = tile.xb;
yaa = tile.ya;
ybb = tile.yb;
tx = ((xbb - xaa) * abscissa + xaa + xbb) / 2;
x = xtrans(tx);
ty = (abscissa' * (ybb - yaa) + yaa + ybb) / 2;
y = ones (15, 1) * ya(x) + ytrans(ty) * (yb(x) - ya(x));
xhalfwidth = (xbb - xaa ) / 2;
yhalfwidth = ones (15, 1) * (yb(x) - ya(x)) .* (ybb - yaa) ./ 2;
x = ones (15, 1) * x;
tx = ones (15,1) * tx;
ty = ty * ones (1, 15);
z = yhalfwidth .* f (x, y) .* ztrans(tx, ty) .* xhalfwidth;
q = weights15 * (weights15 * z)';
qerr = abs(weights7 * (weights7 * z)' - q);
endfunction
## Test input validation
%!error quad2d
%!error quad2d (0, 1 ,2 ,3 ,4)
%!error quad2d (@plus)
%!error quad2d (@plus, 1)
%!error quad2d (@plus, 1, 2)
%!error quad2d (@plus, 1, 2, 3)
%!error quad2d (@plus, 1, 2, 3, [4 5])
%!error quad2d (@plus, 1, 2, 3, "test")
%!error quad2d (@plus, 1, 2, 3, 4, "foo")
%!error quad2d (@plus, 1, 2, 3, 4, "foo", "bar")
%!error quad2d (@plus, 1, 2, 3, 4, 99, "bar")
%!error quad2d (@plus, 1, 2, 3, 4, "AbsTol", "foo")
%!error quad2d (@plus, 1, 2, 3, 4, "AbsTol", [1, 2])
%!error quad2d (@plus, 1, 2, 3, 4, "AbsTol", -1)
%!error quad2d (@plus, 1, 2, 3, 4, "RelTol", "foo")
%!error quad2d (@plus, 1, 2, 3, 4, "RelTol", [1, 2])
%!error quad2d (@plus, 1, 2, 3, 4, "RelTol", -1)
%!warning quad2d (@plus, 1, 2, 3, 4, "MaxFunEvals", 1);
%!test
%! f = @(x, y) x .* y;
%! assert (quad2d (f, 0, 1, 0, 1), 0.25, 1e-10);
%!test
%! f = @(x, y) 9 * x.^2 + 15 * y.^2;
%! assert (quad2d (f, 0, 5, -5, 0, "AbsTol", 1e-9), 5000, 1e-9);
%! assert (quad2d (f, 0, 5, -5, 0, "RelTol", 1e-6), 5000, -1e-6);
%! assert (quad2d (f, 0, 5, -5, 0, "RelTol", 1e-6, "AbsTol", 1e-9),
%! 5000, 1e-9);
## tests from dblquad
%!assert (quad2d (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-7),
%! 2*log (2), 1e-6)
%!assert (quad2d (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "RelTol", 1e-6),
%! 2*log (2), -1e-6)
%!assert (quad2d (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-8,
%! "RelTol", 1e-6), 2*log (2), -1e-6)
%!assert (quad2d (@(x, y) exp (-x.^2 - y.^2) , -1, 1, -1, 1),
%! pi * erf (1).^2, 1e-10)
%!assert (quad2d (@plus, 1, 2, 3, 4), 5, 1e-10)
%!assert (quad2d (@(x,y) 1 ./ (x + y), 0, 1, 0, @(x) 1 - x), 1, -1e-6)
%!assert (quad2d (@(x,y) 1 ./ (x + y), 0, 1, 0, @(x) 1 - x, "Singular", true),
%! 1, -1e-6)