1######################################################################## 2## 3## Copyright (C) 2017-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {@var{q} =} integral3 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb}) 28## @deftypefnx {} {@var{q} =} integral3 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb}, @var{prop}, @var{val}, @dots{}) 29## 30## Numerically evaluate the three-dimensional integral of @var{f} using 31## adaptive quadrature over the three-dimensional domain defined by 32## @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb} (scalars may 33## be finite or infinite). Additionally, @var{ya} and @var{yb} may be 34## scalar functions of @var{x} and @var{za}, and @var{zb} maybe be scalar 35## functions of @var{x} and @var{y}, allowing for integration over 36## non-rectangular domains. 37## 38## @var{f} is a function handle, inline function, or string containing the name 39## of the function to evaluate. The function @var{f} must be of the form 40## @math{z = f(x,y)} where @var{x} is a vector and @var{y} is a scalar. It 41## should return a vector of the same length and orientation as @var{x}. 42## 43## Additional optional parameters can be specified using 44## @qcode{"@var{property}", @var{value}} pairs. Valid properties are: 45## 46## @table @code 47## @item AbsTol 48## Define the absolute error tolerance for the quadrature. The default 49## value is 1e-10 (1e-5 for single). 50## 51## @item RelTol 52## Define the relative error tolerance for the quadrature. The default 53## value is 1e-6 (1e-4 for single). 54## 55## @item Method 56## Specify the two-dimensional integration method to be used, with valid 57## options being @qcode{"auto"} (default), @qcode{"tiled"}, or 58## @qcode{"iterated"}. When using @qcode{"auto"}, Octave will choose the 59## @qcode{"tiled"} method unless any of the integration limits are infinite. 60## 61## @item Vectorized 62## Enable or disable vectorized integration. A value of @code{false} forces 63## Octave to use only scalar inputs when calling the integrand, which enables 64## integrands @math{f(x,y)} that have not been vectorized and only accept 65## @var{x} and @var{y} as scalars to be used. The default value is 66## @code{true}. 67## @end table 68## 69## Adaptive quadrature is used to minimize the estimate of error until the 70## following is satisfied: 71## @tex 72## $$error \leq \max \left( AbsTol, RelTol\cdot\vert q\vert \right)$$ 73## @end tex 74## @ifnottex 75## 76## @example 77## @group 78## @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|) 79## @end group 80## @end example 81## 82## @end ifnottex 83## 84## @var{err} is an approximate bound on the error in the integral 85## @code{abs (@var{q} - @var{I})}, where @var{I} is the exact value of the 86## integral. 87## 88## Example 1 : integrate over a rectangular volume 89## 90## @example 91## @group 92## @var{f} = @@(@var{x},@var{y},@var{z}) ones (size (@var{x})); 93## @var{q} = integral3 (@var{f}, 0, 1, 0, 1, 0, 1) 94## @result{} @var{q} = 1.00000 95## @end group 96## @end example 97## 98## For this constant-value integrand, the result is a volume which is just 99## @code{@var{Length} * @var{Width} * @var{Height}}. 100## 101## Example 2 : integrate over a spherical volume 102## 103## @example 104## @group 105## @var{f} = @@(@var{x},@var{y}) ones (size (@var{x})); 106## @var{ymax} = @@(@var{x}) sqrt (1 - @var{x}.^2); 107## @var{zmax} = @@(@var{x},@var{y}) sqrt (1 - @var{x}.^2 - @var{y}.^2); 108## @var{q} = integral3 (@var{f}, 0, 1, 0, @var{ymax}, 0, @var{zmax}) 109## @result{} @var{q} = 0.52360 110## @end group 111## @end example 112## 113## For this constant-value integrand, the result is a volume which is 1/8th 114## of a unit sphere or @code{1/8 * 4/3 * pi}. 115## 116## Programming Notes: If there are singularities within the integration region 117## it is best to split the integral and place the singularities on the 118## boundary. 119## 120## Known @sc{matlab} incompatibility: If tolerances are left unspecified, and 121## any integration limits are of type @code{single}, then Octave's integral 122## functions automatically reduce the default absolute and relative error 123## tolerances as specified above. If tighter tolerances are desired they 124## must be specified. @sc{matlab} leaves the tighter tolerances appropriate 125## for @code{double} inputs in place regardless of the class of the 126## integration limits. 127## 128## Reference: @nospell{L.F. Shampine}, 129## @cite{@sc{matlab} program for quadrature in 2D}, Applied Mathematics and 130## Computation, pp.@: 266--274, Vol 1, 2008. 131## 132## @seealso{triplequad, integral, quad, quadgk, quadv, quadl, 133## quadcc, trapz, integral2, quad2d, dblquad} 134## @end deftypefn 135 136function q = integral3 (f, xa, xb, ya, yb, za, zb, varargin) 137 138 if (nargin < 7 || mod (nargin, 2) == 0) 139 print_usage (); 140 endif 141 142 if (! is_function_handle (f)) 143 print_usage (); 144 endif 145 146 if (! (isreal (xa) && isscalar (xa) && isreal (xb) && isscalar (xb))) 147 print_usage (); 148 endif 149 150 ## Check for single or double limits to set appropriate default tolerance. 151 issingle = (isa ([xa, xb], "single") 152 || (! is_function_handle (ya) && isa (ya, "single")) 153 || (! is_function_handle (yb) && isa (yb, "single")) 154 || (! is_function_handle (za) && isa (za, "single")) 155 || (! is_function_handle (zb) && isa (zb, "single"))); 156 157 ## Communicate to downstream quadrature routines that at least one limit of 158 ## integration was of single type by casting xa, xb to single. 159 if (issingle) 160 xa = single (xa); 161 xb = single (xb); 162 endif 163 164 ## Set default tolerances, and then update with any specified parameters. 165 if (issingle) 166 abstol = 1e-5; 167 reltol = 1e-4; 168 else 169 abstol = 1e-10; 170 reltol = 1e-6; 171 endif 172 173 method = "auto"; 174 vectorized = true; 175 idx = 1; 176 while (idx < nargin - 7) 177 prop = varargin{idx++}; 178 if (! ischar (prop)) 179 error ("integral3: property PROP must be a string"); 180 endif 181 182 switch (tolower (prop)) 183 case "abstol" 184 abstol = varargin{idx++}; 185 if (! (isnumeric (abstol) && isscalar (abstol) && abstol >= 0)) 186 error ("integral3: AbsTol value must be a numeric scalar >= 0"); 187 endif 188 189 case "reltol" 190 reltol = varargin{idx++}; 191 if (! (isnumeric (reltol) && isscalar (reltol) && reltol >= 0)) 192 error ("integral3: RelTol value must be a numeric scalar >= 0"); 193 endif 194 195 case "method" 196 method = tolower (varargin{idx++}); 197 if (! any (strcmp (method, {"auto", "iterated", "tiled"}))) 198 error ("integral3 : unrecognized method '%s'", method); 199 endif 200 201 case "vectorized" 202 vectorized = varargin{idx++}; 203 if (! (isscalar (vectorized) && isreal (vectorized))) 204 error ('integral3: Vectorized must be a logical value'); 205 endif 206 207 otherwise 208 error ("integral3: unknown property '%s'", prop); 209 210 endswitch 211 endwhile 212 213 if (strcmp (method, "auto")) 214 if (isinf (xa) || isinf (xb) 215 || (! is_function_handle (ya) && isinf (ya)) 216 || (! is_function_handle (yb) && isinf (yb)) 217 || (! is_function_handle (za) && isinf (za)) 218 || (! is_function_handle (zb) && isinf (zb))) 219 method = "iterated"; 220 else 221 method = "tiled"; 222 endif 223 endif 224 225 ## check upper and lower bounds of y 226 if (! is_function_handle (ya)) 227 if (! (isreal (ya) && isscalar (ya))) 228 error ("integral3: YA must be a real scalar or a function"); 229 endif 230 ya = @(x) ya * ones (size (x)); 231 endif 232 if (! is_function_handle (yb)) 233 if (! (isreal (yb) && isscalar (yb))) 234 error ("integral3: YB must be a real scalar or a function"); 235 endif 236 yb = @(x) yb * ones (size (x)); 237 endif 238 239 ## check upper and lower bounds of z 240 if (! is_function_handle (za)) 241 if (! (isreal (za) && isscalar (za))) 242 error ("integral3: ZA must be a real scalar or a function"); 243 endif 244 za = @(x, y) za * ones (size(y)); 245 endif 246 if (! is_function_handle (zb)) 247 if (! (isreal (zb) && isscalar (zb))) 248 error ("integral3: ZB must be a real scalar or a function") 249 endif 250 zb = @(x, y) zb * ones (size (y)); 251 endif 252 253 finner = @(x) inner (x, f, ya, yb, za, zb, vectorized, method, abstol, reltol); 254 q = quadcc (finner, xa, xb, [abstol, reltol]); 255 256endfunction 257 258function q = inner (x, f, ya, yb, za, zb, vectorized, method, abstol, reltol) 259 q = zeros (size (x)); 260 for i = 1 : length (x) 261 za2 = @(y) za(x(i), y); 262 zb2 = @(y) zb(x(i), y); 263 f2 = @(y, z) f(x(i), y, z); 264 if (! vectorized) 265 f2 = @(y, z) arrayfun (f2, y, z); 266 endif 267 if (strcmp (method, "iterated")) 268 finner_iter = @(y) inner_iterated (y, f2, za2, zb2, abstol, reltol); 269 q(i) = quadcc (finner_iter, ya(x(i)), yb(x(i)), [abstol, reltol]); 270 else 271 q(i) = quad2d (f2, ya(x(i)), yb(x(i)), za2, zb2, 272 "AbsTol", abstol, "RelTol", reltol); 273 endif 274 endfor 275endfunction 276 277function q = inner_iterated (y, f2, za2, zb2, abstol, reltol) 278 q = zeros (size (y)); 279 for i = 1 : length (y) 280 q(i) = quadcc (@(z) f2(y(i), z), za2(y(i)), zb2(y(i)), [abstol, reltol]); 281 endfor 282endfunction 283 284 285## method tests 286%!shared f 287%! f = @(x, y, z) x .* y .* z; 288 289%!assert (integral3 (f, 0, 1, 0, 1, 0, 1), 0.125, 1e-10); 290%!assert (integral3 (f, 0, 1, 0, 1, 0, 1, "method", "tiled"), 0.125, 1e-10); 291%!assert (integral3 (f, 0, 1, 0, 1, 0, 1, "method", "iterated"), 0.125, 1e-10); 292%!assert (integral3 (f, 0, 1, 0, 1, 0, 1, "method", "auto"), 0.125, 1e-10); 293 294## vectorized = false test 295%!test 296%! f = @(x, y, z) x * y * z; 297%! assert (integral3 (f, 0, 1, 0, 1, 0, 1, "vectorized", false), 0.125, 1e-10); 298 299## tolerance tests 300%!test 301%! f = @(x, y, z) 2 * x.^2 + 3 * y.^2 + 4 * z.^2; 302%!assert (integral3 (f, 0, 5, -5, 0, 0, 5, "AbsTol", 1e-9), 9375, 1e-9) 303%!assert (integral3 (f, 0, 5, -5, 0, 0, 5, "RelTol", 1e-5), 9375, -1e-5) 304%!assert (integral3 (f, 0, 5, -5, 0, 0, 5, "RelTol", 1e-6, "AbsTol", 1e-9), 305%! 9375, 1e-9) 306 307## non-rectangular region 308## This test is too slow with "iterated" method 309%!test 310%! f = @(x,y,z) 1 ./ (x + y + z); 311%! ymax = @(x) 1 - x; 312%! zmax = @(x, y) 1 - x - y; 313%! assert (integral3 (f, 0, 1, 0, ymax, 0, zmax, "method", "tiled"), 0.25, 1e-6); 314 315## Test input validation 316%!error integral3 317%!error integral3 (@plus) 318%!error integral3 (@plus, 1) 319%!error integral3 (@plus, 1, 2) 320%!error integral3 (@plus, 1, 2, 3) 321%!error integral3 (@plus, 1, 2, 3, 4) 322%!error integral3 (@plus, 1, 2, 3, 4, 5) 323%!error integral3 (@plus, 1, 2, 3, 4, 5, 6, "foo") 324%!error integral3 (0, 1, 2, 3, 4, 5, 6) # f must be a function handle 325%!error integral3 (@plus, 1i, 2, 3, 4, 5, 6) # real limits 326%!error integral3 (@plus, 1, 2i, 3, 4, 5, 6) # real limits 327%!error integral3 (@plus, [1 1], 2, 3, 4, 5, 6) # scalar limits 328%!error integral3 (@plus, 1, [2 2], 3, 4, 5, 6) # scalar limits 329%!error <property PROP must be a string> 330%! integral3 (@plus, 1, 2, 3, 4, 5, 6, 99, "bar"); 331%!error <AbsTol value must be a numeric> 332%! integral3 (@plus, 1, 2, 3, 4, 5, 6, "AbsTol", "foo"); 333%!error <AbsTol value must be a .* scalar> 334%! integral3 (@plus, 1, 2, 3, 4, 5, 6, "AbsTol", [1, 2]); 335%!error <AbsTol value must be.* .= 0> 336%! integral3 (@plus, 1, 2, 3, 4, 5, 6, "AbsTol", -1); 337%!error <RelTol value must be a numeric> 338%! integral3 (@plus, 1, 2, 3, 4, 5, 6, "RelTol", "foo"); 339%!error <RelTol value must be a .* scalar> 340%! integral3 (@plus, 1, 2, 3, 4, 5, 6, "RelTol", [1, 2]); 341%!error <RelTol value must be.* .= 0> 342%! integral3 (@plus, 1, 2, 3, 4, 5, 6, "RelTol", -1); 343%!error <unrecognized method 'foo'> 344%! integral3 (@plus,1,2,3,4,5,6, "method", "foo"); 345%!error <Vectorized must be a logical value> 346%! integral3 (@plus,1,2,3,4,5,6, "Vectorized", [0 1]); 347%!error <Vectorized must be a logical value> 348%! integral3 (@plus,1,2,3,4,5,6, "Vectorized", {true}); 349%!error <unknown property 'foo'> 350%! integral3 (@plus, 1, 2, 3, 4, 6, 6, "foo", "bar"); 351%!error <YA must be a real scalar> integral3 (@plus, 1, 2, 3i, 4, 5, 6) 352%!error <YA must be a real scalar> integral3 (@plus, 1, 2, [3 3], 4, 5, 6) 353%!error <YB must be a real scalar> integral3 (@plus, 1, 2, 3, 4i, 5, 6) 354%!error <YB must be a real scalar> integral3 (@plus, 1, 2, 3, [4 4], 5, 6) 355%!error <ZA must be a real scalar> integral3 (@plus, 1, 2, 3, 4, 5i, 6) 356%!error <ZA must be a real scalar> integral3 (@plus, 1, 2, 3, 4, [5 5], 6) 357%!error <ZB must be a real scalar> integral3 (@plus, 1, 2, 3, 4, 5, 6i) 358%!error <ZB must be a real scalar> integral3 (@plus, 1, 2, 3, 4, 5, [6 6]) 359