## 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} =} integral3 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb}) ## @deftypefnx {} {@var{q} =} integral3 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb}, @var{prop}, @var{val}, @dots{}) ## ## Numerically evaluate the three dimensional integral of @var{f} using adaptive ## quadrature over the three-dimensional domain defined by @var{xa}, @var{xb}, ## @var{ya}, @var{yb}, @var{za}, @var{zb} (scalars may be finite or infinite). ## ## @code{integral3} is a wrapper for @code{triplequad} 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{w = f(x,y,z)} where either @var{x} or @var{y} is a vector and the ## remaining inputs are scalars. @var{f} should return a vector of the same ## length and orientation as the vector input @var{x} or @var{y}. ## ## 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{integral3} currently only functions over a rectangular volume. ## Implementing @var{ya} and @var{yb} as functions of @var{x}, and @var{za} and ## @var{zb} as functions of (@var{x,y})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 3d 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 = integral3 (f, xa, xb, ya, yb, za, zb, 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 ## integral3 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 < 7 || (mod (nargin, 2) == 0)) print_usage (); endif if (! is_function_handle (f)) print_usage (); endif if ((! isscalar (xa)) || (! isscalar (xb)) ... || (! isscalar (ya)) || (! isscalar (yb)) ... || (! isscalar (za)) || (! isscalar (zb))) 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") || ... isa (za, "single") || isa (zb, "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 - 7) prop = varargin{idx++}; if (! ischar (prop)) error ("integral3: property PROP must be a string"); endif switch (tolower (prop)) case "reltol" reltol = varargin{idx++}; if (! ischar (reltol)) warning("integral3: 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("integral3: AbsTol value must be a scalar") endif case "method" intmethod = varargin{idx++}; warning (["integral3: alternate integration methods not yet ", ... "implemented. Method property ignored."]); otherwise error ("integral3: 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 ("integral3: unknown RelTol value '%s'", reltol); endif endif q = triplequad (f, xa, xb, ya, yb, za, zb, abstol, integfunc); endfunction %!test %! f = @(x, y, z) x.*y.*z; %! assert (integral3 (f, 0, 1, 0, 1, 0, 1), 0.125, 1e-10); %!test %! f = @(x,y,z) y.*sin(x)+z.*cos(x); %! assert (integral3 (f, 0, pi, 0, 1, -1, 1), 2, 1e-10); %! assert (integral3 (f, 0, pi, 0, 1, -1, 1, 'RelTol', 'off'), 2, 1e-10); ## tests from triplequad %! assert (integral3 (@(x,y,z) exp (-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1, %! 1), pi^(3/2) * erf (1).^3, 1e-10); %! assert (integral3 (@(x,y,z) exp (-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1, %! 1, "RelTol", "off"), pi^(3/2) * erf (1).^3, 1e-10); ## Test input validation %!error integral3 (0, 1 ,2 ,3 ,4, 5, 6) %!error integral3 (@plus) %!error integral3 (@plus, 1) %!error integral3 (@plus, 1, 2) %!error integral3 (@plus, 1, 2, 3) %!error integral3 (@plus, 1, 2, 3, 4) %!error integral3 (@plus, 1, 2, 3, 4, 5) %!error integral3 (@plus, 1, 2, 3, 4, 5, [6 7]) %!error integral3 (@plus, 1, 2, 3, 4, 5, 'test') %!error integral3 (@plus, 1, 2, 3, 4, 5, 6, 'foo') %!error integral3 (@plus, 1, 2, 3, 4, 5, 6, 'foo', 'bar') %!error integral3 (@plus, 1, 2, 3, 4, 5, 6, 99, 'bar') %!error integral3 (@plus, 1, 2, 3, 4, 5, 6, 'AbsTol', 'foo') %!error integral3 (@plus, 1, 2, 3, 4, 5, 6, 'RelTol', 'foo')