## 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 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. ## ## @item Vectorized ## Option to enable/disable vectorized integration. False forces octave to use ## only scalar inputs when calling the integrand, which enables integrands ## @math{f(x,y)} that have not been vectorized and only accept @var{x} and ## @var{y} as scalars. Default value is @code{true}. ## ## @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 ## ## @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 ## 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, 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 method = "auto"; 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" method = tolower (varargin{idx++}); if (! any (strcmp (method, {"auto", "iterated", "tiled"}))) error ("integral2 : method '%s' unrecognized", method) endif case "vectorized" # option to allow unvectorized functions to be used vectorized = varargin{idx++}; if (! islogical (vectorized)) error ("integral2: 'vectorized' must be a logical value"); elseif (! vectorized) f = @(x, y) arrayfun(f, x, y); endif otherwise error ("integral2: unknown property '%s'", prop); endswitch endwhile if strcmp (method, "auto") if ((isinf (xa)) || (isinf (xb)) || ... ((! is_function_handle(ya)) && (isinf (ya))) || ... ((! is_function_handle(yb)) && (isinf (yb)))) method = "iterated"; else method = "tiled"; endif 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 ("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 if strcmp(method, "iterated") 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 [q, err] = tiled_integral2(f, xa, xb, ya, yb, abstol, reltol); endif endfunction function q = outer_iterated (f, xa, xb, ya, yb, abstol, reltol) 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 function [q, qerr] = tiled_integral2(f, xa, xb, ya, yb, abstol, reltol) maxiter = 5000; iter = 0; qaccept = 0; qerraccept = 0; area = (xb - xa); # Initialize tile list tilelist(1) = struct("xa", xa, "xb", xb, "ya", 0, "yb", 1, "q", 0, "qerr", Inf); 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 tiles(4) = struct("xa", tile.xa, "xb", tile.xa + tile.xb / 2, "ya", tile.ya, "yb", tile.ya + tile.yb / 2, "q", 0, "qerr", 0); tiles(3) = struct("xa", tile.xa, "xb", tile.xa + tile.xb / 2, "ya", tile.ya + tile.yb / 2, "yb", tile.yb, "q", 0, "qerr", 0); tiles(2) = struct("xa", tile.xa + tile.xb / 2, "xb", tile.xb, "ya", tile.ya, "yb", tile.ya + tile.yb / 2, "q", 0, "qerr", 0); tiles(1) = struct("xa", tile.xa + tile.xb / 2, "xb", tile.xb, "ya", tile.ya + tile.yb / 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)); endfor q = qaccept + sum([[tilelist.q], tiles.q]); qerr = qerraccept + sum([[tilelist.qerr], tiles.qerr]); tol = max (abstol, reltol .* abs (q)); localtol = tol * ([tile.xb] - [tile.xa]) * ([tile.yb] - [tile.ya]) / area; # If global tolerance is met, return if (qerr < 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 (qerr > max (abstol, reltol .* abs (q))) warning ("integral2 : Maximum number of sub-tiles without convergence.") else warning("integral2: Maximum number of sub-tiles, maybe low accuracy.") endif endif endfunction function [q, qerr] = tensorproduct (f, ya, yb, tile) # The Shampine TwoD paper proposes using a G3,K7 rule in a tensor product # I couldn't find a tabulated abscissas et weighrs of a G3,K7 rule publically # available, so I'm using a G7,G15 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; xcenter = (xbb + xaa) / 2; xhalfwidth = (xbb - xaa) / 2; x = xhalfwidth * abscissa + xcenter; xx = ones(15, 1) * x; ylower = ya(x); yupper = yb(x); yhalfwidth = (yupper - ylower) .* (ybb - yaa) ./ 2; ycenter = ylower + (yupper - ylower) .* (ybb + yaa) ./ 2; yy = (abscissa' * yhalfwidth + ones(15, 1) * ycenter); zz = f (xx, yy); q = weights15 * (weights15 * (ones(15,1) * yhalfwidth .* zz))' .* xhalfwidth; qerr = abs(weights7 * (weights7 * (ones(15,1) * yhalfwidth .* zz))' .* xhalfwidth - q); endfunction ## method tests %!test %! f = @(x, y) x .* y; %! assert (integral2 (f, 0, 1, 0, 1), 0.25, 1e-10); %! assert (integral2 (f, 0, 1, 0, 1, "method", "tiled"), 0.25, 1e-10); %! assert (integral2 (f, 0, 1, 0, 1, "method", "iterated"), 0.25, 1e-10); %! assert (integral2 (f, 0, 1, 0, 1, "method", "auto"), 0.25, 1e-10); ## vectorized = false test s %!test %! f = @(x, y) x * y; %! assert (integral2 (f, 0, 1, 0, 1, "vectorized", false), 0.25, 1e-10); %!error integral2 (f, 0, 1, 0, 1, "vectorized", true) ## tolerance tests %!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-6); %!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) %!error integral2 (@plus, 1, 2, 3, 4, "method", "bad")