## 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')