Mercurial > octave
changeset 24148:2e64bed0bb3a
Updated integral2.m and integral3.m and added quad2d (bug #52074).
* /scripts/general/quad2d.m: New function file implementing tiled integration.
* integral2.m: Updated to avoid dblquad and make use of quad2d.
* integral3.m: Updated to avoid double or triplequad and make use of quad2d.
* /scripts/general/module.mk: Add quad2d.m to build system.
* /scripts/help/__unimplemented__.m: Removed quad2d from unimplemented list.
Changed BIST test which had relied on quad2d being unimplemented.
* quad.txi: Updated 2D&3D integration section, added quad2d DOCSTRING.
* NEWS: Announce new function quad2d.
| author | Nicholas R. Jankowski <jankowskin@asme.org> |
|---|---|
| date | Sat, 14 Oct 2017 18:30:51 -0400 |
| parents | 9fab3273ff26 |
| children | 1a3229a2f1ab |
| files | NEWS doc/interpreter/quad.txi scripts/general/integral2.m scripts/general/integral3.m scripts/general/module.mk scripts/general/quad2d.m scripts/help/__unimplemented__.m |
| diffstat | 7 files changed, 975 insertions(+), 249 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Sat Oct 14 18:16:55 2017 -0500 +++ b/NEWS Sat Oct 14 18:30:51 2017 -0400 @@ -78,6 +78,7 @@ integral2 integral3 openvar + quad2d repelem rticks thetaticks
--- a/doc/interpreter/quad.txi Sat Oct 14 18:16:55 2017 -0500 +++ b/doc/interpreter/quad.txi Sat Oct 14 18:30:51 2017 -0400 @@ -274,13 +274,12 @@ @node Functions of Multiple Variables @section Functions of Multiple Variables -Octave does not have built-in functions for computing the integral of -functions of multiple variables directly. It is possible, however, to -compute the integral of a function of multiple variables using the -existing functions for one-dimensional integrals. +Octave includes several functions for computing the integral of functions of +multiple variables. This procedure can generally be performed by creating a +function that integrates @math{f} with respect to @math{x}, and then integrates +that function with respect to @math{y}. This procedure can be performed +manually using the following example which integrates the function: -To illustrate how the integration can be performed, we will integrate -the function @tex $$ f(x, y) = \sin(\pi x y)\sqrt{x y} @@ -289,19 +288,16 @@ @ifnottex @example -f(x, y) = sin(pi*x*y)*sqrt(x*y) +f(x, y) = sin(pi*x*y) * sqrt(x*y) @end example @end ifnottex for @math{x} and @math{y} between 0 and 1. -The first approach creates a function that integrates @math{f} with -respect to @math{x}, and then integrates that function with respect to -@math{y}. Because @code{quad} is written in Fortran it cannot be called -recursively. This means that @code{quad} cannot integrate a function -that calls @code{quad}, and hence cannot be used to perform the double -integration. Any of the other integrators, however, can be used which is -what the following code demonstrates. +Using @code{quadgk} in the example below, a double integration can be +performed. (Note that any of the 1-D quadrature functions can be used in this +fashion except for @code{quad} since it is written in Fortran and cannot be +called recursively.) @example @group @@ -318,9 +314,10 @@ @end group @end example -The above process can be simplified with the @code{dblquad} and -@code{triplequad} functions for integrals over two and three variables. For -example: +The algorithm above is implemented in the function @code{dblquad} for integrals +over two variables. The 3-D equivalent of this process is implemented in +@code{triplequad} for integrals over three variables. As an example, the +result above can be replicated with a call to @code{dblquad} as shown below. @example @group @@ -333,19 +330,32 @@ @DOCSTRING(triplequad) -@code{integral2} and @code{integral3} are also provided as partial -compatibility wrappers to @code{dblquad} and @code{triplequad}: +The recursive algorithm for quadrature presented above is referred to as +@qcode{"iterated"}. A separate 2-D integration method is implemented in the +function @code{quad2d}. This function performs a @qcode{"tiled"} integration +by subdividing the integration domain into rectangular regions and performing +separate integrations over those domains. The domains are further subdivided +in areas requiring refinement to reach the desired numerical accuracy. For +certain functions this method can be faster than the 2-D iteration used in the +other functions above. + +@DOCSTRING(quad2d) + +Finally, the functions @code{integral2} and @code{integral3} are provided +as general 2-D and 3-D integration functions. They will auto-select between +iterated and tiled integration methods and, unlike @code{dblquad} and +@code{triplequad}, will work with non-rectangular integration domains. @DOCSTRING(integral2) @DOCSTRING(integral3) -The above mentioned approach works, but is fairly slow, and that problem -increases exponentially with the dimensionality of the integral. Another -possible solution is to use Orthogonal Collocation as described in the -previous section (@pxref{Orthogonal Collocation}). The integral of a function -@math{f(x,y)} for @math{x} and @math{y} between 0 and 1 can be approximated -using @math{n} points by +The above integrations can be fairly slow, and that problem increases +exponentially with the dimensionality of the integral. Another possible +solution for 2-D integration is to use Orthogonal Collocation as described in +the previous section (@pxref{Orthogonal Collocation}). The integral of a +function @math{f(x,y)} for @math{x} and @math{y} between 0 and 1 can be +approximated using @math{n} points by @tex $$ \int_0^1 \int_0^1 f(x,y) d x d y \approx \sum_{i=1}^n \sum_{j=1}^n q_i q_j f(r_i, r_j),
--- a/scripts/general/integral2.m Sat Oct 14 18:16:55 2017 -0500 +++ b/scripts/general/integral2.m Sat Oct 14 18:30:51 2017 -0400 @@ -1,4 +1,4 @@ -## Copyright (C) 2017 Nicholas Jankowski +## Copyright (C) 2017 Nicholas Jankowski, David Bateman ## ## This file is part of Octave. ## @@ -19,14 +19,13 @@ ## -*- 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). -## -## @code{integral2} is a wrapper for @code{dblquad} intended to provide -## @sc{matlab} compatibility. More control of the numerical integration may -## be achievable by calling the various quadrature functions directly. +## @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 +## 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 @@ -44,6 +43,19 @@ ## @item RelTol ## Define the relative error tolerance for the quadrature. The default ## value is 1e-6 (1e-4 for single). +## +## @item Method +## Specify the two-dimensional integration method to be used, with valid +## options being @qcode{"auto"} (default), @qcode{"tiled"}, or +## @qcode{"iterated"}. When using @qcode{"auto"}, Octave will choose the +## @qcode{"tiled"} method unless any of the integration limits are infinite. +## +## @item Vectorized +## Enable or disable vectorized integration. A value of @code{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 to be used. The default value is +## @code{true}. ## @end table ## ## Adaptive quadrature is used to minimize the estimate of error until the @@ -55,59 +67,67 @@ ## ## @example ## @group -## @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|). +## @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|) ## @end group ## @end example ## ## @end ifnottex ## -## Known @sc{matlab} incompatibilities: +## @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. +## +## Example 1 : integrate a rectangular region in x-y plane +## +## @example +## @group +## @var{f} = @@(@var{x},@var{y}) 2*ones (size (@var{x})); +## @var{q} = integral2 (@var{f}, 0, 1, 0, 1) +## @result{} @var{q} = 2 +## @end group +## @end example +## +## The result is a volume, which for this constant-value integrand, is just +## @code{@var{Length} * @var{Width} x @var{Height}}. +## +## Example 2 : integrate a triangular region in x-y plane ## -## @enumerate -## @item -## If tolerances are left unspecified, and any integration limits or waypoints -## 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. -# -## @item -## @code{integral2} currently only integrates functions over rectangular -## domains. Implementing @var{ya} and @var{yb} as functions of @var{x} is a -## planned future improvement. +## @example +## @group +## @var{f} = @@(@var{x},@var{y}) 2*ones (size (@var{x})); +## @var{ymax} = @@(@var{x}) 1 - @var{x}; +## @var{q} = integral2 (@var{f}, 0, 1, 0, @var{ymax}) +## @result{} @var{q} = 1 +## @end group +## @end example +## +## The result is a volume, which for this constant-value integrand, is the +## Triangle Area x Height or +## @code{1/2 * @var{Base} * @var{Width} * @var{Height}}. ## -## @item -## The @qcode{"Method"} property is not yet implemented in Octave due to the -## lack of a @qcode{"tiled"} integrator implementation. All integrals are -## evaluated using an equivalent of the @qcode{"iterated"} method. -## @end enumerate +## Programming Notes: If there are singularities within the integration region +## it is best to split the integral and place the singularities on the +## boundary. +## +## Known @sc{matlab} incompatibility: 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. ## -## @seealso{integral, integral3, quad, quadgk, quadv, quadl, quadcc, trapz, -## dblquad, triplequad} +## Reference: @nospell{L.F. Shampine}, +## @cite{@sc{matlab} program for quadrature in 2D}, Applied Mathematics and +## Computation, pp. 266--274, Vol 1, 2008. +## +## @seealso{quad2d, dblquad, integral, quad, quadgk, quadv, quadl, quadcc, +## trapz, integral3, 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 that - ## 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. +function [q, err] = integral2 (f, xa, xb, ya, yb, varargin) - ## 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 re-create 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. - - if (nargin < 5 || (mod (nargin, 2) == 0)) + if (nargin < 5 || mod (nargin, 2) == 0) print_usage (); endif @@ -115,12 +135,14 @@ print_usage (); endif - if (! (isscalar (xa) && isscalar (xb) && isscalar (ya)) && isscalar (yb)) + if (! (isreal (xa) && isscalar (xa) && isreal (xb) && isscalar (xb))) print_usage (); endif ## Check for single or double limits to set appropriate default tolerance. - issingle = isa ([xa, xb, ya, yb], "single"); + issingle = (isa ([xa, xb], "single") + || (! is_function_handle (ya) && isa (ya, "single")) + || (! is_function_handle (yb) && isa (yb, "single"))); ## Set defaults, update with any specified parameters. if (issingle) @@ -131,91 +153,186 @@ reltol = 1e-6; endif - if (nargin == 5) - ## Pass the simplest case directly to integrator. - ## Let quadcc function handle input checks. - q = dblquad (f, xa, xb, ya, yb, [abstol, reltol], @quadcc); + 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 - else - ## Parse options to determine how to call integrator. - intmethod = []; + case "reltol" + reltol = varargin{idx++}; + if (! (isnumeric (reltol) && isscalar (reltol) && reltol >= 0)) + error ("integral2: RelTol value must be a numeric scalar >= 0"); + endif - idx = 1; - while (idx < nargin - 5) - prop = varargin{idx++}; - if (! ischar (prop)) - error ("integral2: property PROP must be a string"); - endif + case "method" + method = tolower (varargin{idx++}); + if (! any (strcmp (method, {"auto", "iterated", "tiled"}))) + error ("integral2 : unrecognized method '%s'", method); + endif + + case "vectorized" + vectorized = varargin{idx++}; + if (! (isscalar (vectorized) && isreal (vectorized))) + error ('integral2: Vectorized must be a logical value'); + endif + if (! vectorized) + f = @(x, y) arrayfun (f, x, y); + 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 + 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 - case "reltol" - reltol = varargin{idx++}; - if (! (isnumeric (reltol) && isscalar (reltol) && reltol >= 0)) - error ("integral2: RelTol value must be a numeric scalar >= 0"); - endif + ## check upper and lower bounds of y + if (! is_function_handle (ya)) + if (! (isreal (ya) && isscalar (ya))) + error ("integral2: YA must be a real scalar or a function"); + endif + ya = @(x) ya * ones (rows (x), columns (x)); + endif + if (! is_function_handle (yb)) + if (! (isreal (yb) && isscalar (yb))) + error ("integral2: YB must be a real scalar or a function"); + endif + yb = @(x) yb * ones (rows (x), columns (x)); + endif - case "method" - intmethod = varargin{idx++}; - warning (["integral2: Only 'iterated' method implemented. ", ... - "Method property ignored."]); - otherwise - error ("integral2: unknown property '%s'", prop); - endswitch - endwhile + if (strcmp (method, "iterated")) + q = outer_iterated (f, xa, xb, ya, yb, abstol, reltol); - q = dblquad (f, xa, xb, ya, yb, [abstol, reltol], @quadcc); - + if (nargout == 2) + warning ('integral2: "iterated" method can not return estimated error'); + err = 0; + endif + else + [q, err] = quad2d (f, xa, xb, ya, yb, "AbsTol", abstol, "RelTol", reltol); endif endfunction +function q = outer_iterated (f, xa, xb, ya, yb, abstol, reltol) + finner_iter = @(x) inner_iterated (x, f, ya, yb, abstol, reltol); + q = quadcc (finner_iter, xa, xb, [abstol, reltol]); +endfunction -%!test +function q = inner_iterated (x, f, ya, yb, abstol, reltol) + q = zeros (size (x)); + for i = 1 : length (x) + q(i) = quadcc (@(y) f(x(i), y), ya(x(i)), yb(x(i)), [abstol, reltol]); + endfor +endfunction + + +## method tests +%!shared f %! f = @(x, y) x .* y; -%! assert (integral2 (f, 0, 1, 0, 1), 0.25, 1e-10); +%!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 %!test +%! f = @(x, y) x * y; +%! assert (integral2 (f, 0, 1, 0, 1, "vectorized", false), 0.25, 1e-10); + +## tolerance tests +%!shared f %! 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); + +%!assert (integral2 (f, 0, 5, -5, 0, "AbsTol", 1e-9), 5000, 1e-9) +%!assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-5), 5000, -1e-5) +%!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-7) -%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "RelTol", 1e-6), +%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "RelTol", 1e-5), +%! 2*log (2), -1e-5) +%!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) 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) + +## tests from dblquad w/method specified +%!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-5, "method", "iterated"), +%! 2*log (2), -1e-5) +%!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 () %!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 (0, 1, 2, 3, 4) # f must be function handle +%!error integral2 (@plus, 1i, 2, 3, 4) # real limits +%!error integral2 (@plus, 1, 2i, 3, 4) # real limits +%!error integral2 (@plus, [1 1], 2, 3, 4) # scalar limits +%!error integral2 (@plus, 1, [2 2], 3, 4) # scalar limits +%!error <property PROP must be a string> integral2 (@plus,1,2,3,4,99, "bar") +%!error <AbsTol value must be a numeric> +%! integral2 (@plus,1,2,3,4, "AbsTol", "foo"); +%!error <AbsTol value must be a .* scalar> +%! integral2 (@plus, 1, 2, 3, 4, "AbsTol", [1, 2]); +%!error <AbsTol value must be.* .= 0> integral2 (@plus,1,2,3,4, "AbsTol", -1) +%!error <RelTol value must be a numeric> +%! integral2 (@plus, 1, 2, 3, 4, "RelTol", "foo"); +%!error <RelTol value must be a .* scalar> +%! integral2 (@plus, 1, 2, 3, 4, "RelTol", [1, 2]); +%!error <RelTol value must be.* .= 0> integral2 (@plus,1,2,3,4, "RelTol", -1) +%!error <unrecognized method 'foo'> integral2 (@plus,1,2,3,4, "method", "foo") +%!error <Vectorized must be a logical value> +%! integral2 (@plus,1,2,3,4, "Vectorized", [0 1]); +%!error <Vectorized must be a logical value> +%! integral2 (@plus,1,2,3,4, "Vectorized", {true}); %!error <unknown property 'foo'> integral2 (@plus, 1, 2, 3, 4, "foo", "bar") -%!error <property PROP must be a string> integral2 (@plus, 1, 2, 3, 4, 99, "bar") -%!error <AbsTol value must be a numeric> integral2 (@plus, 1, 2, 3, 4, "AbsTol", "foo") -%!error <AbsTol value must be a .* scalar> integral2 (@plus, 1, 2, 3, 4, "AbsTol", [1, 2]) -%!error <AbsTol value must be.* .= 0> integral2 (@plus, 1, 2, 3, 4, "AbsTol", -1) -%!error <RelTol value must be a numeric> integral2 (@plus, 1, 2, 3, 4, "RelTol", "foo") -%!error <RelTol value must be a .* scalar> integral2 (@plus, 1, 2, 3, 4, "RelTol", [1, 2]) -%!error <RelTol value must be.* .= 0> integral2 (@plus, 1, 2, 3, 4, "RelTol", -1) -%!warning <Only 'iterated' method implemented> -%! q = integral2 (@plus, 0, 1, 0, 1, "Method", "tiled"); +%!error <YA must be a real scalar> integral2 (@plus, 1, 2, 3i, 4) +%!error <YA must be a real scalar> integral2 (@plus, 1, 2, [3 3], 4) +%!error <YB must be a real scalar> integral2 (@plus, 1, 2, 3, 4i) +%!error <YB must be a real scalar> integral2 (@plus, 1, 2, 3, [4 4]) +%!warning <"iterated" method can not return estimated error> +%! [q, err] = integral2 (@plus, 0, 0, 0, 0, "method", "iterated");
--- a/scripts/general/integral3.m Sat Oct 14 18:16:55 2017 -0500 +++ b/scripts/general/integral3.m Sat Oct 14 18:30:51 2017 -0400 @@ -1,4 +1,4 @@ -## Copyright (C) 2017 Nicholas Jankowski +## Copyright (C) 2017 David Bateman ## ## This file is part of Octave. ## @@ -23,17 +23,15 @@ ## 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 -## @sc{matlab} compatibility. More control of the numerical integration may be -## achievable by calling the various quadrature functions directly. +## be finite or infinite). Additionally, @var{ya} and @var{yb} may be +## scalar functions of @var{x} and @var{za}, and @var{zb} maybe be scalar +## functions of @var{x} and @var{y}, allowing for 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{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}. +## @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: @@ -46,6 +44,19 @@ ## @item RelTol ## Define the relative error tolerance for the quadrature. The default ## value is 1e-6 (1e-4 for single). +## +## @item Method +## Specify the two-dimensional integration method to be used, with valid +## options being @qcode{"auto"} (default), @qcode{"tiled"}, or +## @qcode{"iterated"}. When using @qcode{"auto"}, Octave will choose the +## @qcode{"tiled"} method unless any of the integration limits are infinite. +## +## @item Vectorized +## Enable or disable vectorized integration. A value of @code{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 to be used. The default value is +## @code{true}. ## @end table ## ## Adaptive quadrature is used to minimize the estimate of error until the @@ -57,61 +68,67 @@ ## ## @example ## @group -## @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|). +## @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|) ## @end group ## @end example ## ## @end ifnottex ## -## Known @sc{matlab} incompatibilities: +## @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. +## +## Example 1 : integrate over a rectangular volume ## -## @enumerate -## @item -## If tolerances are left unspecified, and any integration limits or waypoints -## 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. +## @example +## @group +## @var{f} = @@(@var{x},@var{y},@var{z}) ones (size (@var{x})); +## @var{q} = integral3 (@var{f}, 0, 1, 0, 1, 0, 1) +## @result{} @var{q} = 1 +## @end group +## @end example +## +## For this constant-value integrand, the result is a volume which is just +## @code{@var{Length} * @var{Width} x @var{Height}}. +## +## Example 2 : integrate over a spherical volume ## -## @item -## @code{integral3} currently only integrates functions over rectangular -## volumes. 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. +## @example +## @group +## @var{f} = @@(@var{x},@var{y}) ones (size (@var{x})); +## @var{ymax} = @@(@var{x}) sqrt (1 - @var{x}.^2); +## @var{zmax} = @@(@var{x}) sqrt (1 - @var{x}.^2 - @var{y}.^2); +## @var{q} = integral3 (@var{f}, 0, 1, 0, @var{ymax}) +## @result{} @var{q} = 0.52360 +## @end group +## @end example +## +## For this constant-value integrand, the result is a volume which is 1/8th +## of a unit sphere or @code{1/8 * 4/3 * pi}. ## -## @item -## The @qcode{"Method"} property is not yet implemented in Octave due to the -## lack of a @qcode{"tiled"} integrator implementation. All integrals are -## evaluated using an equivalent of the @qcode{"iterated"} method. -## @end enumerate +## Programming Notes: If there are singularities within the integration region +## it is best to split the integral and place the singularities on the +## boundary. ## -## @seealso{integral, integral2, quad, quadgk, quadv, quadl, quadcc, trapz, -## dblquad, triplequad} +## Known @sc{matlab} incompatibility: 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. +## +## Reference: @nospell{L.F. Shampine}, +## @cite{@sc{matlab} program for quadrature in 2D}, Applied Mathematics and +## Computation, pp. 266--274, Vol 1, 2008. +## +## @seealso{triplequad, integral, quad, quadgk, quadv, quadl, +## quadcc, trapz, integral2, quad2d, dblquad} ## @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 that - ## 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 underlying 2D integration method to be used on the y and z - ## dimensions, with valid options being @var{"auto"}, @var{"tiled"}, or - ## @var{"iterated"}. @code{integral3} will use @var{"auto"} by default, where - ## it will usually choose @var{"tiled"} unless any of the integration limits - ## are infinite. - - if (nargin < 7 || (mod (nargin, 2) == 0)) + if (nargin < 7 || mod (nargin, 2) == 0) print_usage (); endif @@ -119,14 +136,16 @@ print_usage (); endif - if (! (isscalar (xa) && isscalar (xb) - && isscalar (ya) && isscalar (yb) - && isscalar (za) && isscalar (zb))) + if (! (isreal (xa) && isscalar (xa) && isreal (xb) && isscalar (xb))) print_usage (); endif ## Check for single or double limits to set appropriate default tolerance. - issingle = isa ([xa, xb, ya, yb, za, zb], "single"); + issingle = (isa ([xa, xb], "single") + || (! is_function_handle (ya) && isa (ya, "single")) + || (! is_function_handle (yb) && isa (yb, "single")) + || (! is_function_handle (za) && isa (za, "single")) + || (! is_function_handle (zb) && isa (zb, "single"))); ## Set defaults, update with any specified parameters. if (issingle) @@ -137,81 +156,187 @@ reltol = 1e-6; endif - if (nargin == 7) - ## Pass the simplest case directly to integrator. - ## Let quadcc function handle input checks. - q = triplequad (f, xa, xb, ya, yb, za, zb, [abstol, reltol], @quadcc); + method = "auto"; + vectorized = true; + 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 "abstol" + abstol = varargin{idx++}; + if (! (isnumeric (abstol) && isscalar (abstol) && abstol >= 0)) + error ("integral3: AbsTol value must be a numeric scalar >= 0"); + endif - else - ## Parse options to determine how to call integrator. - intmethod = []; + case "reltol" + reltol = varargin{idx++}; + if (! (isnumeric (reltol) && isscalar (reltol) && reltol >= 0)) + error ("integral3: RelTol value must be a numeric scalar >= 0"); + endif + + case "method" + method = tolower (varargin{idx++}); + if (! any (strcmp (method, {"auto", "iterated", "tiled"}))) + error ("integral3 : unrecognized method '%s'", method); + endif - idx = 1; - while (idx < nargin - 7) - prop = varargin{idx++}; - if (! ischar (prop)) - error ("integral3: property PROP must be a string"); - endif + case "vectorized" + vectorized = varargin{idx++}; + if (! (isscalar (vectorized) && isreal (vectorized))) + error ('integral3: Vectorized must be a logical value'); + endif + + otherwise + error ("integral3: unknown property '%s'", prop); + + endswitch + endwhile - switch (tolower (prop)) - case "abstol" - abstol = varargin{idx++}; - if (! (isnumeric (abstol) && isscalar (abstol) && abstol >= 0)) - error ("integral3: 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 + if (strcmp (method, "auto")) + if (isinf (xa) || isinf (xb) + || (! is_function_handle (ya) && isinf (ya)) + || (! is_function_handle (yb) && isinf (yb)) + || (! is_function_handle (za) && isinf (za)) + || (! is_function_handle (zb) && isinf (zb))) + method = "iterated"; + else + method = "tiled"; + endif + endif - case "method" - intmethod = varargin{idx++}; - warning (["integral3: Only 'iterated' method implemented. ", ... - "Method property ignored."]); - otherwise - error ("integral3: unknown property '%s'", prop); - endswitch - endwhile + ## check upper and lower bounds of y + if (! is_function_handle (ya)) + if (! (isreal (ya) && isscalar (ya))) + error ("integral3: YA must be a real scalar or a function"); + endif + ya = @(x) ya * ones (size (x)); + endif + if (! is_function_handle (yb)) + if (! (isreal (yb) && isscalar (yb))) + error ("integral3: YB must be a real scalar or a function"); + endif + yb = @(x) yb * ones (size (x)); + endif - q = triplequad (f, xa, xb, ya, yb, za, zb, [abstol, reltol], @quadcc); + ## check upper and lower bounds of z + if (! is_function_handle (za)) + if (! (isreal (za) && isscalar (za))) + error ("integral3: ZA must be a real scalar or a function"); + endif + za = @(x, y) za * ones (size(y)); + endif + if (! is_function_handle (zb)) + if (! (isreal (zb) && isscalar (zb))) + error ("integral3: ZB must be a real scalar or a function") + endif + zb = @(x, y) zb * ones (size (y)); + endif - endif + finner = @(x) inner (x, f, ya, yb, za, zb, vectorized, method, abstol, reltol); + q = quadcc (finner, xa, xb, [abstol, reltol]); endfunction +function q = inner (x, f, ya, yb, za, zb, vectorized, method, abstol, reltol) + q = zeros (size (x)); + for i = 1 : length (x) + za2 = @(y) za(x(i), y); + zb2 = @(y) zb(x(i), y); + f2 = @(y, z) f(x(i), y, z); + if (! vectorized) + f2 = @(y, z) arrayfun (f2, y, z); + endif + if (strcmp (method, "iterated")) + finner_iter = @(y) inner_iterated (y, f2, za2, zb2, abstol, reltol); + q(i) = quadcc (finner_iter, ya(x(i)), yb(x(i)), [abstol, reltol]); + else + q(i) = quad2d (f2, ya(x(i)), yb(x(i)), za2, zb2, + "AbsTol", abstol, "RelTol", reltol); + endif + endfor +endfunction -%!test -%! f = @(x, y, z) x.*y.*z; -%! assert (integral3 (f, 0, 1, 0, 1, 0, 1), 0.125, 1e-10); +function q = inner_iterated (y, f2, za2, zb2, abstol, reltol) + q = zeros (size (y)); + for i = 1 : length (y) + q(i) = quadcc (@(z) f2(y(i), z), za2(y(i)), zb2(y(i)), [abstol, reltol]); + endfor +endfunction + +## method tests +%!shared f +%! f = @(x, y, z) x .* y .* z; + +%!assert (integral3 (f, 0, 1, 0, 1, 0, 1), 0.125, 1e-10); +%!assert (integral3 (f, 0, 1, 0, 1, 0, 1, "method", "tiled"), 0.125, 1e-10); +%!assert (integral3 (f, 0, 1, 0, 1, 0, 1, "method", "iterated"), 0.125, 1e-10); +%!assert (integral3 (f, 0, 1, 0, 1, 0, 1, "method", "auto"), 0.125, 1e-10); + +## vectorized = false test %!test -%! f = @(x,y,z) y.*sin(x) + z.*cos(x); -%! assert (integral3 (f, 0, pi, 0, 1, -1, 1), 2, 1e-10); +%! f = @(x, y, z) x * y * z; +%! assert (integral3 (f, 0, 1, 0, 1, 0, 1, "vectorized", false), 0.125, 1e-10); + +## tolerance tests +%!shared f +%! f = @(x, y, z) 2 * x.^2 + 3 * y.^2 + 4 * z.^2; -## 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 (f, 0, 5, -5, 0, 0, 5, "AbsTol", 1e-9), 9375, 1e-9) +%!assert (integral3 (f, 0, 5, -5, 0, 0, 5, "RelTol", 1e-5), 9375, -1e-5) +%!assert (integral3 (f, 0, 5, -5, 0, 0, 5, "RelTol", 1e-6, "AbsTol", 1e-9), +%! 9375, 1e-9) + +## non-rectangular region +## This test is too slow with "iterated" method +%!assert (integral3 (@(x,y,z) 1 ./ (x + y + z), 0, 1, 0, @(x) 1 - x, 0, +%! @(x, y) 1 - x - y, "method", "tiled"), 0.25, 1e-6) ## Test input validation -%!error integral3 (0, 1 ,2 ,3 ,4, 5, 6) +%!error integral3 %!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 <unknown property 'foo'> integral3 (@plus, 1, 2, 3, 4, 5, 6, "foo", "bar") -%!error <property PROP must be a string> integral3 (@plus, 1, 2, 3, 4, 5, 6, NA, "bar") -%!error <AbsTol value must be a numeric> integral3 (@plus, 1, 2, 3, 4, 5, 6, "AbsTol", "foo") -%!error <AbsTol value must be a .* scalar> integral3 (@plus, 1, 2, 3, 4, 5, 6, "AbsTol", [1, 2]) -%!error <AbsTol value must be.* .= 0> integral3 (@plus, 1, 2, 3, 4, 5, 6, "AbsTol", -1) -%!error <RelTol value must be a numeric> integral3 (@plus, 1, 2, 3, 4, 5, 6, "RelTol", "foo") -%!error <RelTol value must be a .* scalar> integral3 (@plus, 1, 2, 3, 4, 5, 6, "RelTol", [1, 2]) -%!error <RelTol value must be.* .= 0> integral3 (@plus, 1, 2, 3, 4, 5, 6, "RelTol", -1) -%!warning <Only 'iterated' method implemented> -%! q = integral3 (@plus, 0, 1, 0, 1, 0, 1, "Method", "tiled"); +%!error integral3 (0, 1, 2, 3, 4, 5, 6) # f must be a function handle +%!error integral3 (@plus, 1i, 2, 3, 4, 5, 6) # real limits +%!error integral3 (@plus, 1, 2i, 3, 4, 5, 6) # real limits +%!error integral3 (@plus, [1 1], 2, 3, 4, 5, 6) # scalar limits +%!error integral3 (@plus, 1, [2 2], 3, 4, 5, 6) # scalar limits +%!error <property PROP must be a string> +%! integral3 (@plus, 1, 2, 3, 4, 5, 6, 99, "bar"); +%!error <AbsTol value must be a numeric> +%! integral3 (@plus, 1, 2, 3, 4, 5, 6, "AbsTol", "foo"); +%!error <AbsTol value must be a .* scalar> +%! integral3 (@plus, 1, 2, 3, 4, 5, 6, "AbsTol", [1, 2]); +%!error <AbsTol value must be.* .= 0> +%! integral3 (@plus, 1, 2, 3, 4, 5, 6, "AbsTol", -1); +%!error <RelTol value must be a numeric> +%! integral3 (@plus, 1, 2, 3, 4, 5, 6, "RelTol", "foo"); +%!error <RelTol value must be a .* scalar> +%! integral3 (@plus, 1, 2, 3, 4, 5, 6, "RelTol", [1, 2]); +%!error <RelTol value must be.* .= 0> +%! integral3 (@plus, 1, 2, 3, 4, 5, 6, "RelTol", -1); +%!error <unrecognized method 'foo'> +%! integral3 (@plus,1,2,3,4,5,6, "method", "foo"); +%!error <Vectorized must be a logical value> +%! integral3 (@plus,1,2,3,4,5,6, "Vectorized", [0 1]); +%!error <Vectorized must be a logical value> +%! integral3 (@plus,1,2,3,4,5,6, "Vectorized", {true}); +%!error <unknown property 'foo'> +%! integral3 (@plus, 1, 2, 3, 4, 6, 6, "foo", "bar"); +%!error <YA must be a real scalar> integral3 (@plus, 1, 2, 3i, 4, 5, 6) +%!error <YA must be a real scalar> integral3 (@plus, 1, 2, [3 3], 4, 5, 6) +%!error <YB must be a real scalar> integral3 (@plus, 1, 2, 3, 4i, 5, 6) +%!error <YB must be a real scalar> integral3 (@plus, 1, 2, 3, [4 4], 5, 6) +%!error <ZA must be a real scalar> integral3 (@plus, 1, 2, 3, 4, 5i, 6) +%!error <ZA must be a real scalar> integral3 (@plus, 1, 2, 3, 4, [5 5], 6) +%!error <ZB must be a real scalar> integral3 (@plus, 1, 2, 3, 4, 5, 6i) +%!error <ZB must be a real scalar> integral3 (@plus, 1, 2, 3, 4, 5, [6 6])
--- a/scripts/general/module.mk Sat Oct 14 18:16:55 2017 -0500 +++ b/scripts/general/module.mk Sat Oct 14 18:30:51 2017 -0400 @@ -66,6 +66,7 @@ %reldir%/postpad.m \ %reldir%/prepad.m \ %reldir%/publish.m \ + %reldir%/quad2d.m \ %reldir%/quadgk.m \ %reldir%/quadl.m \ %reldir%/quadv.m \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/general/quad2d.m Sat Oct 14 18:30:51 2017 -0400 @@ -0,0 +1,476 @@ +## 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 +## <http://www.gnu.org/licenses/>. + +## -*- 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 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). +## +## @item MaxFunEvals +## The maximum number of function calls to the vectorized function @var{f}. +## The default value is 5000. +## +## @item Singular +## Enable/disable transforms to weaken singularities on the edge of the +## integration domain. 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 FailurePlot +## If @code{quad2d} fails to converge to the desired error tolerance before +## MaxFunEvals is reached, a plot of the areas that still need refinement +## is created. The default value is @var{false}. +## @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. The optional output @var{iter} is the number of vectorized +## function calls to the function @var{f} that were used. +## +## Example 1 : integrate a rectangular region in x-y plane +## +## @example +## @group +## @var{f} = @@(@var{x},@var{y}) 2*ones (size (@var{x})); +## @var{q} = quad2d (@var{f}, 0, 1, 0, 1) +## @result{} @var{q} = 2 +## @end group +## @end example +## +## The result is a volume, which for this constant-value integrand, is just +## @code{@var{Length} * @var{Width} x @var{Height}}. +## +## Example 2 : integrate a triangular region in x-y plane +## +## @example +## @group +## @var{f} = @@(@var{x},@var{y}) 2*ones (size (@var{x})); +## @var{ymax} = @@(@var{x}) 1 - @var{x}; +## @var{q} = quad2d (@var{f}, 0, 1, 0, @var{ymax}) +## @result{} @var{q} = 1 +## @end group +## @end example +## +## The result is a volume, which for this constant-value integrand, is the +## Triangle Area x Height or +## @code{1/2 * @var{Base} * @var{Width} * @var{Height}}. +## +## Programming Notes: If there are singularities within the integration region +## it is best to split the integral and place the singularities on the +## boundary. +## +## Known @sc{matlab} incompatibility: 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. +## +## Reference: @nospell{L.F. Shampine}, +## @cite{@sc{matlab} program for quadrature in 2D}, Applied Mathematics and +## Computation, pp. 266--274, Vol 1, 2008. +## +## @seealso{integral2, dblquad, integral, quad, quadgk, quadv, quadl, quadcc, +## trapz, integral3, triplequad} +## @end deftypefn + +function [q, err, iter] = quad2d (f, xa, xb, ya, yb, varargin) + + if (nargin < 5 || mod (nargin, 2) == 0) + print_usage (); + endif + + if (ischar (f)) + ## Convert function given as a string to a function handle + f = @(x) feval (f, x); + elseif (! is_function_handle (f)) + print_usage (); + endif + + if (! (isreal (xa) && isscalar (xa) && isreal (xb) && isscalar (xb))) + print_usage (); + endif + + ## Check for single or double limits to set appropriate default tolerance. + issingle = (isa ([xa, xb], "single") + || (! is_function_handle (ya) && isa (ya, "single")) + || (! 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 "maxfunevals" + maxiter = varargin{idx++}; + if (! (isnumeric (maxiter) && isscalar (maxiter) + && fix (maxiter) == maxiter && maxiter >= 1)) + error ("quad2d: MaxFunEvals value must be a scalar integer >= 1"); + endif + + case "singular" + singular = varargin{idx++}; + if (! (isscalar (singular) && isreal (singular))) + error ("quad2d: Singular must be a logical value"); + endif + + case "vectorized" + vectorized = varargin{idx++}; + if (! (isscalar (vectorized) && isreal (vectorized))) + error ("quad2d: Vectorized must be a logical value"); + endif + + case "failureplot" + failureplot = varargin{idx++}; + if (! (isscalar (failureplot) && isreal (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 (! (isreal (ya) && isscalar (ya))) + error ("quad2d: YA must be a real scalar or a function"); + endif + ya = @(x) ya * ones(rows (x), columns (x)); + endif + if (! is_function_handle (yb)) + if (! (isreal (yb) && isscalar (yb))) + error ("quad2d: YB must be a real scalar or a function"); + endif + yb = @(x) yb * ones(rows (x), columns (x)); + 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 + [~, idx] = max ([tilelist.qerr]); + tile = tilelist(idx); + tilelist(idx) = []; + + ## 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. + idx = find ([tiles.qerr] < localtol); + qaccept += sum ([tiles(idx).q]); + qerraccept += sum ([tiles(idx).qerr]); + tiles(idx) = []; + tilelist = [tilelist, tiles]; + endwhile + + ## Verify convergence + if (iter >= maxiter) + if (err > max (abstol, reltol .* abs (q))) + warning (["quad2d: " ... + "Maximum number of sub-tiles reached 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 publicly + ## available, so use 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 + + +%!shared f +%! 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 +%!test +%! f = @(x, y) 1 ./ (x+y); +%!assert (quad2d (f, 0, 1, 0, 1, "AbsTol", 1e-7), 2*log (2), 1e-7) +%!assert (quad2d (f, 0, 1, 0, 1, "RelTol", 1e-5), 2*log (2), -1e-5) +%!assert (quad2d (f, 0, 1, 0, 1, "AbsTol", 1e-8, "RelTol", 1e-6), +%! 2*log (2), -1e-6) +%!assert (quad2d (f, 0, 1, 0, @(x) 1 - x), 1, -1e-6) +%!assert (quad2d (f, 0, 1, 0, @(x) 1 - x, "Singular", true), 1, -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) + +## Test input validation +%!error quad2d () +%!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, "foo") +%!error quad2d (0, 1, 2, 3, 4) # f must be function handle +%!error quad2d (@plus, 1i, 2, 3, 4) # real limits +%!error quad2d (@plus, 1, 2i, 3, 4) # real limits +%!error quad2d (@plus, [1 1], 2, 3, 4) # scalar limits +%!error quad2d (@plus, 1, [2 2], 3, 4) # scalar limits +%!error <property PROP must be a string> quad2d (@plus, 1, 2, 3, 4, 99, "bar") +%!error <AbsTol value must be a numeric> quad2d (@plus, 1, 2, 3, 4, "AbsTol", "foo") +%!error <AbsTol value must be a .* scalar> quad2d (@plus, 1, 2, 3, 4, "AbsTol", [1, 2]) +%!error <AbsTol value must be.* .= 0> quad2d (@plus, 1, 2, 3, 4, "AbsTol", -1) +%!error <RelTol value must be a numeric> quad2d (@plus, 1, 2, 3, 4, "RelTol", "foo") +%!error <RelTol value must be a .* scalar> quad2d (@plus, 1, 2, 3, 4, "RelTol", [1, 2]) +%!error <RelTol value must be.* .= 0> quad2d (@plus, 1, 2, 3, 4, "RelTol", -1) +%!error <MaxFunEvals value must be a scalar integer> +%! quad2d (@plus,1,2,3,4, "MaxFunEvals", {1}); +%!error <MaxFunEvals value must be a scalar integer> +%! quad2d (@plus,1,2,3,4, "MaxFunEvals", [1 1]); +%!error <MaxFunEvals value must be a scalar integer> +%! quad2d (@plus,1,2,3,4, "MaxFunEvals", 1.5); +%!error <MaxFunEvals value must be a scalar integer .= 1> +%! quad2d (@plus,1,2,3,4, "MaxFunEvals", -1); +%!error <Singular must be a logical value> +%! quad2d (@plus,1,2,3,4, "Singular", [0 1]); +%!error <Singular must be a logical value> +%! quad2d (@plus,1,2,3,4, "Singular", {true}); +%!error <Vectorized must be a logical value> +%! quad2d (@plus,1,2,3,4, "Vectorized", [0 1]); +%!error <Vectorized must be a logical value> +%! quad2d (@plus,1,2,3,4, "Vectorized", {true}); +%!error <FailurePlot must be a logical value> +%! quad2d (@plus,1,2,3,4, "FailurePlot", [0 1]); +%!error <FailurePlot must be a logical value> +%! quad2d (@plus,1,2,3,4, "FailurePlot", {true}); +%!error <unknown property 'foo'> quad2d (@plus, 1, 2, 3, 4, "foo", "bar") +%!error <YA must be a real scalar> quad2d (@plus, 1, 2, 3i, 4) +%!error <YA must be a real scalar> quad2d (@plus, 1, 2, [3 3], 4) +%!error <YB must be a real scalar> quad2d (@plus, 1, 2, 3, 4i) +%!error <YB must be a real scalar> quad2d (@plus, 1, 2, 3, [4 4]) +%!warning <Maximum number of> quad2d (@plus, 1, 2, 3, 4, "MaxFunEvals", 1);
--- a/scripts/help/__unimplemented__.m Sat Oct 14 18:16:55 2017 -0500 +++ b/scripts/help/__unimplemented__.m Sat Oct 14 18:30:51 2017 -0400 @@ -66,7 +66,7 @@ case "matlabrc" txt = ["matlabrc is not implemented. ", ... - 'Octave uses the file ".octaverc" instead.']; + "Octave uses the file '.octaverc' instead."]; case {"ode113", "ode15i", "ode15s", "ode23s", "ode23t", "ode23tb"} txt = ["Octave provides lsode and ode45 for solving differential equations. ", ... @@ -78,9 +78,6 @@ txt = ["startup is not implemented. ", ... 'Octave uses the file ".octaverc" instead.']; - case "quad2d" - txt = ["quad2d is not implemented. Consider using dblquad."]; - case {"xlsread", "xlsfinfo", "xlswrite", "wk1read", "wk1finfo", "wk1write"} txt = ["Functions for spreadsheet style I/O ", ... "(.xls .xlsx .sxc .ods .dbf .wk1 etc.) " , ... @@ -843,11 +840,10 @@ }; endfunction - %!test %! str = __unimplemented__ ("no_name_function"); %! assert (isempty (str)); -%! str = __unimplemented__ ("quad2d"); -%! assert (str(1:51), "quad2d is not implemented. Consider using dblquad."); +%! str = __unimplemented__ ("matlabrc"); +%! assert (str(1:71), "matlabrc is not implemented. Octave uses the file '.octaverc' instead."); %! str = __unimplemented__ ("MException"); %! assert (str(1:58), "the 'MException' function is not yet implemented in Octave");