Mercurial > octave

comparison scripts/general/quadgk.m @ 31058:12f8fb75fc30

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
quadgk.m: Overhaul function to add new "ArrayValued" option (bug #62468) * quadgk.m: Add "ArrayValued" option to docstring. Re-phrase bits of existing documentation for clarity. Rename input validation variable "str" to "prop" for clarity. Add decode for "arrayvalued" property. If limits of integration are high to low then reverse them in function rather than recursively calling quadgk(). Add FIXME note about missing input validation for several properties. Delete variable 'h0' which was never used. New code branch if "ArrayValued" is true which uses same code strategy as for normal quadgk, but subfunction __quadgk_eval_array__ for evaluating the function in a vectorized manner. Add BIST tests for "ArrayValued" input. * quadgk.m (__quadgk_eval__): Eliminate "too_close" output and code to calculate "too_close" which was never used. * quadgk.m (__quadgk_eval_array__): New subfunction.
author Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org>
date Thu, 02 Jun 2022 08:45:10 -0700
parents 4298af839d20
children 451fb63a10a0
comparison
equal deleted inserted replaced
31057:b437597ec652 31058:12f8fb75fc30
25 25
26 ## -*- texinfo -*- 26 ## -*- texinfo -*-
27 ## @deftypefn {} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}) 27 ## @deftypefn {} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b})
28 ## @deftypefnx {} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}, @var{abstol}) 28 ## @deftypefnx {} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}, @var{abstol})
29 ## @deftypefnx {} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}, @var{abstol}, @var{trace}) 29 ## @deftypefnx {} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}, @var{abstol}, @var{trace})
30 ## @deftypefnx {} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}, @var{prop}, @var{val}, @dots{}) 30 ## @deftypefnx {} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}, "@var{prop}", @var{val}, @dots{})
31 ## @deftypefnx {} {[@var{q}, @var{err}] =} quadgk (@dots{}) 31 ## @deftypefnx {} {[@var{q}, @var{err}] =} quadgk (@dots{})
32 ## 32 ##
33 ## Numerically evaluate the integral of @var{f} from @var{a} to @var{b} 33 ## Numerically evaluate the integral of @var{f} from @var{a} to @var{b}
34 ## using adaptive @nospell{Gauss-Kronrod} quadrature. 34 ## using adaptive @nospell{Gauss-Kronrod} quadrature.
35 ## 35 ##
36 ## @var{f} is a function handle, inline function, or string containing the name 36 ## @var{f} is a function handle, inline function, or string containing the name
37 ## of the function to evaluate. The function @var{f} must be vectorized and 37 ## of the function to evaluate. The function @var{f} must be vectorized and
38 ## return a vector of output values when given a vector of input values. 38 ## return a vector of output values when given a vector of input values (See
39 ## property @qcode{"ArrayValued"} for an exception to this rule).
39 ## 40 ##
40 ## @var{a} and @var{b} are the lower and upper limits of integration. Either 41 ## @var{a} and @var{b} are the lower and upper limits of integration. Either
41 ## or both limits may be infinite or contain weak end singularities. Variable 42 ## or both limits may be infinite or contain weak end singularities. Variable
42 ## transformation will be used to treat any infinite intervals and weaken the 43 ## transformation will be used to treat any infinite intervals and weaken the
43 ## singularities. For example: 44 ## singularities. For example:
58 ## interval and evaluating each subinterval. If @var{trace} is true then after 59 ## interval and evaluating each subinterval. If @var{trace} is true then after
59 ## computing each of these partial integrals display: (1) the number of 60 ## computing each of these partial integrals display: (1) the number of
60 ## subintervals at this step, (2) the current estimate of the error @var{err}, 61 ## subintervals at this step, (2) the current estimate of the error @var{err},
61 ## (3) the current estimate for the integral @var{q}. 62 ## (3) the current estimate for the integral @var{q}.
62 ## 63 ##
63 ## The behavior of the algorithm can be configured by passing arguments 64 ## The behavior of the algorithm can be configured by passing arguments to
64 ## to @code{quadgk} as pairs @qcode{"@var{prop}", @var{val}}. Valid properties 65 ## @code{quadgk} as pairs @qcode{"@var{prop}", @var{val}}. Valid properties
65 ## are 66 ## are
66 ## 67 ##
67 ## @table @code 68 ## @table @code
68 ## @item AbsTol 69 ## @item AbsTol
69 ## Define the absolute error tolerance for the quadrature. The default 70 ## Define the absolute error tolerance for the quadrature. The default
71 ## 72 ##
72 ## @item RelTol 73 ## @item RelTol
73 ## Define the relative error tolerance for the quadrature. The default 74 ## Define the relative error tolerance for the quadrature. The default
74 ## relative tolerance is 1e-6 (1e-4 for single). 75 ## relative tolerance is 1e-6 (1e-4 for single).
75 ## 76 ##
77 ## @item ArrayValued
78 ## When set to true, the function @var{f} produces an array output for a scalar
79 ## input. The default is false which requires that @var{f} produce an output
80 ## that is the same size as the input. For example,
81 ##
82 ## @example
83 ## quadgk (@@(x) x .^ (1:5), 0, 2, "ArrayValued", 1)
84 ## @end example
85 ##
86 ## will integrate @code{[x.^1, x.^2, x.^3, x.^4, x.^5]} in one function call
87 ## rather than having to repeatedly define a single anonymous function and
88 ## use a normal invocation of @code{quadgk}.
89 ##
90 ## @item WayPoints
91 ## Specify points which will become endpoints for subintervals in the
92 ## algorithm which can result in significantly improved estimation of the error
93 ## in the integral, faster computation, or both. It can be useful to specify
94 ## more subintervals around a region where the integrand is rapidly changing or
95 ## to flag locations where there is a discontinuity in the first derivative
96 ## of the function. For example, the signum function has a discontinuity at
97 ## @code{x == 0} and by specifying a waypoint
98 ##
99 ## @example
100 ## quadgk (@@(x) sign (x), -0.5, 1, "Waypoints", [0])
101 ## @end example
102 ##
103 ## @noindent
104 ## the error bound is reduced from 4e-7 to 1e-13.
105 ##
106 ## If the function has @strong{singularities} within the region of integration
107 ## those should not be addressed with waypoints. Instead, the overall integral
108 ## should be decomposed into a sum of several smaller integrals such that the
109 ## singularity occurs as one of the bounds of integration in the call to
110 ## @code{quadgk}.
111 ##
112 ## If any of the waypoints are complex then contour integration is performed as
113 ## documented below.
114 ##
76 ## @item MaxIntervalCount 115 ## @item MaxIntervalCount
77 ## @code{quadgk} initially subdivides the interval on which to perform the 116 ## @code{quadgk} initially subdivides the interval on which to perform the
78 ## quadrature into 10 intervals. Subintervals that have an unacceptable error 117 ## quadrature into 10 intervals or, if WayPoints are given, at each waypoint.
79 ## are subdivided and re-evaluated. If the number of subintervals exceeds 650 118 ## Subintervals that have an unacceptable error are subdivided and
80 ## subintervals at any point then a poor convergence is signaled and the 119 ## re-evaluated. If the number of subintervals exceeds 650 subintervals at any
81 ## current estimate of the integral is returned. The property 120 ## point then a poor convergence is signaled and the current estimate of the
82 ## @qcode{"MaxIntervalCount"} can be used to alter the number of subintervals 121 ## integral is returned. The property @qcode{"MaxIntervalCount"} can be used
83 ## that can exist before exiting. 122 ## to alter the number of subintervals that can exist before exiting.
84 ##
85 ## @item WayPoints
86 ## Discontinuities in the first derivative of the function to integrate can be
87 ## flagged with the @qcode{"WayPoints"} property. This forces the ends of a
88 ## subinterval to fall on the breakpoints of the function and can result in
89 ## significantly improved estimation of the error in the integral, faster
90 ## computation, or both. For example,
91 ##
92 ## @example
93 ## quadgk (@@(x) abs (1 - x.^2), 0, 2, "Waypoints", 1)
94 ## @end example
95 ##
96 ## @noindent
97 ## signals the breakpoint in the integrand at @code{@var{x} = 1}.
98 ## 123 ##
99 ## @item Trace 124 ## @item Trace
100 ## If logically true @code{quadgk} prints information on the convergence of the 125 ## If logically true @code{quadgk} prints information on the convergence of the
101 ## quadrature at each iteration. 126 ## quadrature at each iteration.
102 ## @end table 127 ## @end table
103 ## 128 ##
104 ## If any of @var{a}, @var{b}, or @var{waypoints} is complex then the 129 ## If any of @var{a}, @var{b}, or @var{waypoints} is complex then the
105 ## quadrature is treated as a contour integral along a piecewise continuous 130 ## quadrature is treated as a contour integral along a piecewise linear
106 ## path defined by 131 ## path defined by
107 ## @code{[@var{a}, @var{waypoints}(1), @var{waypoints}(2), @dots{}, @var{b}]}. 132 ## @code{[@var{a}, @var{waypoints}(1), @var{waypoints}(2), @dots{}, @var{b}]}.
108 ## In this case the integral is assumed to have no edge singularities. For 133 ## In this case the integral is assumed to have no edge singularities. For
109 ## example, 134 ## example,
110 ## 135 ##
134 ## @cite{"Vectorized adaptive quadrature in @sc{matlab}"}, Journal of 159 ## @cite{"Vectorized adaptive quadrature in @sc{matlab}"}, Journal of
135 ## Computational and Applied Mathematics, pp.@: 131--140, Vol 211, Issue 2, 160 ## Computational and Applied Mathematics, pp.@: 131--140, Vol 211, Issue 2,
136 ## Feb 2008. 161 ## Feb 2008.
137 ## 162 ##
138 ## @seealso{quad, quadv, quadl, quadcc, trapz, dblquad, triplequad, integral, 163 ## @seealso{quad, quadv, quadl, quadcc, trapz, dblquad, triplequad, integral,
139 ## integral2, integral3} 164 ## integral2, integral3}
140 ## @end deftypefn 165 ## @end deftypefn
141 166
142 function [q, err] = quadgk (f, a, b, varargin) 167 function [q, err] = quadgk (f, a, b, varargin)
143 168
144 if (nargin < 3) 169 if (nargin < 3)
145 print_usage (); 170 print_usage ();
146 endif
147
148 if (b < a)
149 ## Reverse integration
150 [q, err] = quadgk (f, b, a, varargin{:});
151 q = -q;
152 return;
153 endif 171 endif
154 172
155 abstol = []; 173 abstol = [];
156 reltol = []; 174 reltol = [];
157 waypoints = []; 175 waypoints = [];
158 maxint = 650; 176 maxint = 650;
177 arrayvalued = false;
159 trace = false; 178 trace = false;
160 179
161 ## Parse options if present. 180 ## Parse options if present.
162 if (nargin > 3) 181 if (nargin > 3)
163 if (! ischar (varargin{1})) 182 if (! ischar (varargin{1}))
179 idx = 1; 198 idx = 1;
180 while (idx < nargin - 3) 199 while (idx < nargin - 3)
181 if (! ischar (varargin{idx})) 200 if (! ischar (varargin{idx}))
182 error ("quadgk: property PROP must be a string"); 201 error ("quadgk: property PROP must be a string");
183 endif 202 endif
184 str = varargin{idx++}; 203 prop = varargin{idx++};
185 switch (tolower (str)) 204 switch (tolower (prop))
186 case "reltol" 205 case "reltol"
187 reltol = varargin{idx++}; 206 reltol = varargin{idx++};
188 case "abstol" 207 case "abstol"
189 abstol = varargin{idx++}; 208 abstol = varargin{idx++};
190 case "waypoints" 209 case "waypoints"
191 waypoints = varargin{idx++}(:); 210 waypoints = varargin{idx++}(:);
192 if (isreal (waypoints))
193 waypoints(waypoints < a | waypoints > b) = [];
194 endif
195 case "maxintervalcount" 211 case "maxintervalcount"
196 maxint = varargin{idx++}; 212 maxint = varargin{idx++};
213 case "arrayvalued"
214 arrayvalued = varargin{idx++};
197 case "trace" 215 case "trace"
198 trace = varargin{idx++}; 216 trace = varargin{idx++};
199 otherwise 217 otherwise
200 error ("quadgk: unknown property '%s'", str); 218 error ("quadgk: unknown property '%s'", prop);
201 endswitch 219 endswitch
202 endwhile 220 endwhile
203 endif 221 endif
222 endif
223
224 reverse = 1;
225 contour = iscomplex (a) || iscomplex (b) || iscomplex (waypoints);
226 if ((b < a) && ! contour)
227 ## Reverse integration
228 [b, a] = deal (a, b);
229 waypoints = sort (waypoints(waypoints > a & waypoints < b));
230 reverse = -1;
204 endif 231 endif
205 232
206 issingle = (isa (a, "single") || isa (b, "single") 233 issingle = (isa (a, "single") || isa (b, "single")
207 || isa (waypoints, "single")); 234 || isa (waypoints, "single"));
208 235
216 reltol = ifelse (issingle, 1e-4, 1e-6); 243 reltol = ifelse (issingle, 1e-4, 1e-6);
217 elseif (! isscalar (reltol) || reltol < 0) 244 elseif (! isscalar (reltol) || reltol < 0)
218 error ("quadgk: RELTOL must be a scalar >=0"); 245 error ("quadgk: RELTOL must be a scalar >=0");
219 endif 246 endif
220 247
248 ## FIXME: No validation of inputs MaxIntervalCount, Waypoints, ArrayValued,
249 ## Trace.
250
221 ## Convert function given as a string to a function handle 251 ## Convert function given as a string to a function handle
222 if (ischar (f)) 252 if (ischar (f))
223 f = @(x) feval (f, x); 253 f = @(x) feval (f, x);
224 endif 254 endif
225 255
226 ## Use variable substitution to weaken endpoint singularities and 256 ## Use variable substitution to weaken endpoint singularities and
227 ## to perform integration with endpoints at infinity. 257 ## to perform integration with endpoints at infinity.
228 ## No transform for contour integrals. 258 ## No transform for contour integrals.
229 if (iscomplex (a) || iscomplex (b) || iscomplex (waypoints)) 259 if (contour)
230 ## contour integral, no transform 260 ## contour integral, no transform
231 subs = [a; waypoints; b]; 261 subs = [a; waypoints; b];
232 h = sum (abs (diff (subs))); 262 h = sum (abs (diff (subs)));
233 h0 = h;
234 trans = @(t) t; 263 trans = @(t) t;
264 ## Ensure f is always vectorized even if specified as, e.g., f = @(x) 1;
265 f = @(t) f (t) + 0*t;
235 elseif (isinf (a) && isinf (b)) 266 elseif (isinf (a) && isinf (b))
236 ## Standard infinite to finite integral transformation. 267 ## Standard infinite to finite integral transformation.
237 ## \int_{-\infinity_^\infinity f(x) dx = \int_-1^1 f (g(t)) g'(t) dt 268 ## \int_{-\infinity_^\infinity f(x) dx = \int_-1^1 f (g(t)) g'(t) dt
238 ## where 269 ## where
239 ## g(t) = t / (1 - t^2) 270 ## g(t) = t / (1 - t^2)
245 subs = [-1; trans(waypoints); 1]; 276 subs = [-1; trans(waypoints); 1];
246 else 277 else
247 subs = linspace (-1, 1, 11)'; 278 subs = linspace (-1, 1, 11)';
248 endif 279 endif
249 h = 2; 280 h = 2;
250 h0 = b - a;
251 trans = @(t) t ./ (1 - t.^2); 281 trans = @(t) t ./ (1 - t.^2);
252 f = @(t) f (t ./ (1 - t .^ 2)) .* (1 + t .^ 2) ./ ((1 - t .^ 2) .^ 2); 282 f = @(t) f (t ./ (1 - t .^ 2)) .* (1 + t .^ 2) ./ ((1 - t .^ 2) .^ 2);
253 elseif (isinf (a)) 283 elseif (isinf (a))
254 ## Formula defined in Shampine paper as two separate steps. 284 ## Formula defined in Shampine paper as two separate steps.
255 ## One to weaken singularity at finite end, then a second to transform to 285 ## One to weaken singularity at finite end, then a second to transform to
271 subs = [-1; trans(tmp); 0]; 301 subs = [-1; trans(tmp); 0];
272 else 302 else
273 subs = linspace (-1, 0, 11)'; 303 subs = linspace (-1, 0, 11)';
274 endif 304 endif
275 h = 1; 305 h = 1;
276 h0 = b - a;
277 trans = @(t) b - (t ./ (1 + t)).^2; 306 trans = @(t) b - (t ./ (1 + t)).^2;
278 f = @(s) - 2 * s .* f (b - (s ./ (1 + s)) .^ 2) ./ ((1 + s) .^ 3); 307 f = @(s) - 2 * s .* f (b - (s ./ (1 + s)) .^ 2) ./ ((1 + s) .^ 3);
279 elseif (isinf (b)) 308 elseif (isinf (b))
280 ## Formula defined in Shampine paper as two separate steps. 309 ## Formula defined in Shampine paper as two separate steps.
281 ## One to weaken singularity at finite end, then a second to transform to 310 ## One to weaken singularity at finite end, then a second to transform to
296 subs = [0; trans(tmp); 1]; 325 subs = [0; trans(tmp); 1];
297 else 326 else
298 subs = linspace (0, 1, 11)'; 327 subs = linspace (0, 1, 11)';
299 endif 328 endif
300 h = 1; 329 h = 1;
301 h0 = b - a;
302 trans = @(t) a + (t ./ (1 - t)).^2; 330 trans = @(t) a + (t ./ (1 - t)).^2;
303 f = @(s) 2 * s .* f (a + (s ./ (1 - s)) .^ 2) ./ ((1 - s) .^ 3); 331 f = @(s) 2 * s .* f (a + (s ./ (1 - s)) .^ 2) ./ ((1 - s) .^ 3);
304 else 332 else
305 ## Davis, Rabinowitz, "Methods of Numerical Integration" p441 2ed. 333 ## Davis, Rabinowitz, "Methods of Numerical Integration" p441 2ed.
306 ## Presented in section 5 of the Shampine paper as 334 ## Presented in section 5 of the Shampine paper as
318 subs = [-1; trans(waypoints, a, b); 1]; 346 subs = [-1; trans(waypoints, a, b); 1];
319 else 347 else
320 subs = linspace (-1, 1, 11)'; 348 subs = linspace (-1, 1, 11)';
321 endif 349 endif
322 h = 2; 350 h = 2;
323 h0 = b - a;
324 trans = @(t) ((b - a) ./ 4) * t .* (3 - t.^2) + (b + a) ./ 2; 351 trans = @(t) ((b - a) ./ 4) * t .* (3 - t.^2) + (b + a) ./ 2;
325 f = @(t) f((b - a) ./ 4 .* t .* (3 - t.^2) + (b + a) ./ 2) .* ... 352 f = @(t) f((b - a) ./ 4 .* t .* (3 - t.^2) + (b + a) ./ 2) .* ...
326 3 .* (b - a) ./ 4 .* (1 - t.^2); 353 3 .* (b - a) ./ 4 .* (1 - t.^2);
327 endif 354 endif
328 355
329 ## Split interval into at least 10 subinterval with a 15 point 356 ## Split interval into at least 10 subinterval with a 15 point
330 ## Gauss-Kronrod rule giving a minimum of 150 function evaluations. 357 ## Gauss-Kronrod rule giving a minimum of 150 function evaluations.
331 while (length (subs) < 11) 358 while (numel (subs) < 11)
332 subs = [subs.' ; subs(1:end-1).' + diff(subs.') ./ 2, NaN](:)(1 : end - 1); 359 subs = [subs.' ; subs(1:end-1).' + diff(subs.') ./ 2, NaN](:)(1:end-1);
333 endwhile 360 endwhile
334 subs = [subs(1:end-1), subs(2:end)]; 361 subs = [subs(1:end-1), subs(2:end)];
335 362
336 warn_id = "Octave:quadgk:warning-termination"; 363 warn_id = "Octave:quadgk:warning-termination";
337 364
338 if (issingle) 365 if (! arrayvalued)
339 eps1 = eps ("single"); 366 ## Initial evaluation of the integrand on the subintervals.
367 [q_subs, q_errs] = __quadgk_eval__ (f, subs, trans);
368 q0 = sum (q_subs);
369 err0 = sum (q_errs);
370
371 first = true;
372 while (true)
373 ## Quit if any evaluations are not finite (Inf or NaN).
374 if (any (! isfinite (q_subs)))
375 warning (warn_id, "quadgk: non-finite integrand encountered");
376 q = q0;
377 err = err0;
378 break;
379 endif
380
381 tol = max (abstol, reltol .* abs (q0));
382
383 ## If the global error estimate is met then exit.
384 if (err0 <= tol)
385 q = q0;
386 err = err0;
387 break;
388 endif
389
390 ## Accept the subintervals that meet the convergence criteria.
391 idx = find (abs (q_errs) < tol .* abs (diff (subs, 1, 2)) ./ h);
392 if (first)
393 q = sum (q_subs(idx));
394 err = sum (q_errs(idx));
395 first = false;
396 else
397 q0 = q + sum (q_subs);
398 err0 = err + sum (q_errs);
399 q += sum (q_subs(idx));
400 err += sum (q_errs(idx));
401 endif
402 subs(idx,:) = [];
403
404 ## If no remaining subintervals then exit.
405 if (isempty (subs))
406 break;
407 endif
408
409 if (trace)
410 disp ([rows(subs), err, q0]);
411 endif
412
413 ## Split remaining subintervals in two
414 mid = (subs(:,1) + subs(:,2)) ./ 2;
415 subs = [subs(:,1), mid; mid, subs(:,2)];
416
417 ## If the maximum subinterval count is met, then
418 ## accept remaining subinterval and exit.
419 if (rows (subs) > maxint)
420 warning (warn_id, "quadgk: maximum interval count (%d) exceeded", maxint);
421 q += sum (q_subs);
422 err += sum (q_errs);
423 break;
424 endif
425
426 ## Evaluation of the integrand on the remaining subintervals
427 [q_subs, q_errs] = __quadgk_eval__ (f, subs, trans);
428 endwhile
429
430 if (nargout < 2 && err > max (abstol, reltol * abs (q)))
431 warning (warn_id,
432 "quadgk: Error tolerance not met. Estimated error %g", err);
433 endif
434
435 ## Reverse integral if necessary.
436 q = reverse * q;
437
340 else 438 else
341 eps1 = eps ("double"); 439 ## f is array-valued
342 endif 440 sz = size (f (subs(1)));
343 441
344 ## Initial evaluation of the integrand on the subintervals 442 ## Initial evaluation of the integrand on the subintervals
345 [q_subs, q_errs] = __quadgk_eval__ (f, subs, eps1, trans); 443 [q_subs, q_errs] = __quadgk_eval_array__ (f, subs, trans, prod (sz));
346 q0 = sum (q_subs); 444 q0 = sum (q_subs, 1);
347 err0 = sum (q_errs); 445 err0 = sum (q_errs, 1);
348 446
349 first = true; 447 first = true;
350 while (true) 448 while (true)
351 ## Quit if any evaluations are not finite (Inf or NaN). 449 ## Quit if any evaluations are not finite (Inf or NaN).
352 if (any (! isfinite (q_subs))) 450 if (any (! isfinite (q_subs)(:)))
353 warning (warn_id, "quadgk: non-finite integrand encountered"); 451 warning (warn_id, "quadgk: non-finite integrand encountered");
354 q = q0; 452 q = q0;
355 err = err0; 453 err = err0;
356 break; 454 break;
455 endif
456
457 tol = max (abstol, reltol .* abs (q0));
458
459 ## If the global error estimate is met then exit
460 if (err0 <= tol)
461 q = q0;
462 err = err0;
463 break;
464 endif
465
466 ## Accept subintervals that meet the convergence criteria in all entries.
467 idx = find (all (abs (q_errs) < tol .* abs (diff (subs, 1, 2)) ./ h, 2));
468 if (first)
469 q = sum (q_subs(idx,:), 1);
470 err = sum (q_errs(idx,:), 1);
471 first = false;
472 else
473 q0 = q + sum (q_subs, 1);
474 err0 = err + sum (q_errs, 1);
475 q += sum (q_subs(idx,:), 1);
476 err += sum (q_errs(idx,:), 1);
477 endif
478 subs(idx,:) = [];
479
480 ## If no remaining subintervals exit
481 if (isempty (subs))
482 break;
483 endif
484
485 if (trace)
486 disp ([rows(subs), err(1, 1), q0(1, 1)]); # print only first entry
487 endif
488
489 ## Split remaining subintervals in two
490 mid = (subs(:,1) + subs(:,2)) ./ 2;
491 subs = [subs(:,1), mid; mid, subs(:,2)];
492
493 ## If the maximum subinterval count is met accept remaining subinterval
494 ## and exit
495 if (rows (subs) > maxint)
496 warning (warn_id, "quadgk: maximum interval count (%d) exceeded", maxint);
497 q += sum (q_subs, 1);
498 err += sum (q_errs, 1);
499 break;
500 endif
501
502 ## Evaluation of the integrand on the remaining subintervals
503 [q_subs, q_errs] = __quadgk_eval_array__ (f, subs, trans, prod (sz));
504 endwhile
505
506 i = find (err > max (abstol, reltol * abs (q)), 1);
507 if (nargout < 2 && length (i) > 0)
508 ## like ind2sub, only as vector.
509 j = mod (floor ((i-1)./cumprod ([1 sz(1:end-1)])),sz)+1;
510 s = ["(" sprintf("%d,",j)(1:end-1) ")"];
511 warning (warn_id,
512 "quadgk: Error tolerance not met. First entry at index %s with estimated error %g", s, err(i));
357 endif 513 endif
358 514
359 tol = max (abstol, reltol .* abs (q0)); 515 q = reverse * reshape (q, sz);
360 516 err = reshape (err, sz);
361 ## If the global error estimate is met then exit
362 if (err0 <= tol)
363 q = q0;
364 err = err0;
365 break;
366 endif
367
368 ## Accept the subintervals that meet the convergence criteria.
369 idx = find (abs (q_errs) < tol .* abs (diff (subs, [], 2)) ./ h);
370 if (first)
371 q = sum (q_subs(idx));
372 err = sum (q_errs(idx));
373 first = false;
374 else
375 q0 = q + sum (q_subs);
376 err0 = err + sum (q_errs);
377 q += sum (q_subs(idx));
378 err += sum (q_errs(idx));
379 endif
380 subs(idx,:) = [];
381
382 ## If no remaining subintervals exit
383 if (rows (subs) == 0)
384 break;
385 endif
386
387 if (trace)
388 disp ([rows(subs), err, q0]);
389 endif
390
391 ## Split remaining subintervals in two
392 mid = (subs(:,2) + subs(:,1)) ./ 2;
393 subs = [subs(:,1), mid; mid, subs(:,2)];
394
395 ## If the maximum subinterval count is met accept remaining subinterval
396 ## and exit
397 if (rows (subs) > maxint)
398 warning (warn_id, "quadgk: maximum interval count (%d) exceeded", maxint);
399 q += sum (q_subs);
400 err += sum (q_errs);
401 break;
402 endif
403
404 ## Evaluation of the integrand on the remaining subintervals
405 [q_subs, q_errs] = __quadgk_eval__ (f, subs, eps1, trans);
406 endwhile
407
408 if (nargout < 2 && err > max (abstol, reltol * abs (q)))
409 warning (warn_id,
410 "quadgk: Error tolerance not met. Estimated error %g", err);
411 endif 517 endif
412 518
413 endfunction 519 endfunction
414 520
415 ## FIXME: too_close output is never used in function that calls this one. 521 function [q, err] = __quadgk_eval__ (f, subs, trans)
416 function [q, err, too_close] = __quadgk_eval__ (f, subs, eps1, trans)
417 522
418 ## A (15,7) point pair of Gauss-Kronrod quadrature rules. 523 ## A (15,7) point pair of Gauss-Kronrod quadrature rules.
419 ## The abscissa and weights are copied directly from dqk15w.f from quadpack. 524 ## The abscissa and weights are copied directly from dqk15w.f from quadpack.
420 525
421 persistent abscissa = [-0.9914553711208126e+00, -0.9491079123427585e+00, ... 526 persistent abscissa = [-0.9914553711208126e+00, -0.9491079123427585e+00, ...
441 diag ([0.1294849661688697e+00, 0.2797053914892767e+00, ... 546 diag ([0.1294849661688697e+00, 0.2797053914892767e+00, ...
442 0.3818300505051889e+00, 0.4179591836734694e+00, ... 547 0.3818300505051889e+00, 0.4179591836734694e+00, ...
443 0.3818300505051889e+00, 0.2797053914892767e+00, ... 548 0.3818300505051889e+00, 0.2797053914892767e+00, ...
444 0.1294849661688697e+00]); 549 0.1294849661688697e+00]);
445 550
446 halfwidth = diff (subs, [], 2) ./ 2; 551 halfwidth = diff (subs, 1, 2) ./ 2;
447 center = sum (subs, 2) ./ 2; 552 center = sum (subs, 2) ./ 2;
448 t = (halfwidth * abscissa) + center; 553 t = (halfwidth * abscissa) + center;
449 x = trans ([t(:,1), t(:,end)]); 554 x = trans ([t(:,1), t(:,end)]);
450 555
451 ## Shampine suggests 100 * eps1, beginning of section 6.
452 if (any (abs (diff (x, [], 2) ./ max (abs (x), [], 2))) < 100 * eps1)
453 too_close = true;
454 q = 0;
455 err = 0;
456 return;
457 endif
458
459 too_close = false;
460 y = reshape (f (t(:)), size (t)); 556 y = reshape (f (t(:)), size (t));
461 557
462 ## This is faster than using bsxfun as the * operator can use a 558 ## This is faster than using bsxfun as the * operator can use a
463 ## single BLAS call, rather than rows (sub) calls to the @times function. 559 ## single BLAS call, rather than rows (sub) calls to the @times function.
464 q = sum (y * weights15, 2) .* halfwidth; 560 q = sum (y * weights15, 2) .* halfwidth;
465 err = abs (sum (y(:,2:2:end) * weights7, 2) .* halfwidth - q); 561 err = abs (sum (y(:,2:2:end) * weights7, 2) .* halfwidth - q);
562
563 endfunction
564
565 function [q, err] = __quadgk_eval_array__ (f, subs, trans, nel)
566
567 ## A (15,7) point pair of Gauss-Kronrod quadrature rules.
568 ## The abscissa and weights are copied directly from dqk15w.f from quadpack.
569
570 persistent abscissa = [-0.9914553711208126e+00, -0.9491079123427585e+00, ...
571 -0.8648644233597691e+00, -0.7415311855993944e+00, ...
572 -0.5860872354676911e+00, -0.4058451513773972e+00, ...
573 -0.2077849550078985e+00, 0.0000000000000000e+00, ...
574 0.2077849550078985e+00, 0.4058451513773972e+00, ...
575 0.5860872354676911e+00, 0.7415311855993944e+00, ...
576 0.8648644233597691e+00, 0.9491079123427585e+00, ...
577 0.9914553711208126e+00];
578
579 persistent weights15 = ...
580 [0.2293532201052922e-01, 0.6309209262997855e-01, ...
581 0.1047900103222502e+00, 0.1406532597155259e+00, ...
582 0.1690047266392679e+00, 0.1903505780647854e+00, ...
583 0.2044329400752989e+00, 0.2094821410847278e+00, ...
584 0.2044329400752989e+00, 0.1903505780647854e+00, ...
585 0.1690047266392679e+00, 0.1406532597155259e+00, ...
586 0.1047900103222502e+00, 0.6309209262997855e-01, ...
587 0.2293532201052922e-01];
588
589 persistent weights7 = ...
590 [0.1294849661688697e+00, 0.2797053914892767e+00, ...
591 0.3818300505051889e+00, 0.4179591836734694e+00, ...
592 0.3818300505051889e+00, 0.2797053914892767e+00, ...
593 0.1294849661688697e+00];
594
595 halfwidth = diff (subs, 1, 2) ./ 2;
596 center = sum (subs, 2) ./ 2;
597 t = (halfwidth * abscissa) + center;
598 x = trans ([t(:,1), t(:,end)]);
599
600 y = zeros (nel, columns(t), rows(t));
601 for i = 1:rows (t)
602 for j = 1:columns(t)
603 y(:,j,i) = f (t(i,j))(:);
604 endfor
605 endfor
606 y = permute (y, [2 3 1]);
607
608 q = reshape (weights15 * y(:,:), [rows(t), nel]) .* halfwidth;
609 err = abs (reshape (weights7 * y(2:2:end,:), rows (t), nel) .* halfwidth - q);
466 610
467 endfunction 611 endfunction
468 612
469 function t = __quadgk_finite_waypoint__ (x, a, b) 613 function t = __quadgk_finite_waypoint__ (x, a, b)
470 c = (-4 .* x + 2.* (b + a)) ./ (b - a); 614 c = (-4 .* x + 2.* (b + a)) ./ (b - a);
500 %!assert (quadgk (@(x) exp (-x .^ 2),-Inf,0), sqrt (pi)/2, -1e-6) 644 %!assert (quadgk (@(x) exp (-x .^ 2),-Inf,0), sqrt (pi)/2, -1e-6)
501 %!test 645 %!test
502 %! f = @(x) x .^ 5 .* exp (-x) .* sin (x); 646 %! f = @(x) x .^ 5 .* exp (-x) .* sin (x);
503 %! assert (quadgk (f, 0, Inf, "RelTol", 1e-8, "AbsTol", 1e-12), -15, -1e-8); 647 %! assert (quadgk (f, 0, Inf, "RelTol", 1e-8, "AbsTol", 1e-12), -15, -1e-8);
504 648
505 %!test <*62412> 649 ## Test vector-valued functions
650 %!assert (quadgk (@(x) [(sin (x)), (sin (2 * x))], 0, pi, "arrayvalued", 1),
651 %! [2, 0], 1e-6)
652
653 ## Test matrix-valued functions
654 %!assert (quadgk (@(x) [ x,x,x; x,1./sqrt(x),x; x,x,x ], 0, 1, "arrayvalued",1),
655 %! [0.5,0.5,0.5; 0.5,2,0.5; 0.5,0.5,0.5], 15*1e-6);
656
657 ## Bug #62412
658 %!warning <Error tolerance not met>
506 %! f = @(t) -1 ./ t.^1.1; 659 %! f = @(t) -1 ./ t.^1.1;
507 %! fail ("quadgk (f, 1, Inf)", "warning", "Error tolerance not met"); 660 %! quadgk (f, 1, Inf);
508 %! [q, err] = quadgk (f, 1, Inf);
509 %! assert (err > 1e-5);
510 661
511 ## Test input validation 662 ## Test input validation
512 %!error quadgk (@sin) 663 %!error quadgk (@sin)
513 %!error quadgk (@sin, 0) 664 %!error quadgk (@sin, 0)
514 %!error <can not pass additional arguments> quadgk (@sin, 0, 1, 1e-6, true, 4) 665 %!error <can not pass additional arguments> quadgk (@sin, 0, 1, 1e-6, true, 4)