Mercurial > octave
view scripts/general/quadv.m @ 29358:0a5b15007766 stable
update Octave Project Developers copyright for the new year
In files that have the "Octave Project Developers" copyright notice,
update for 2021.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 10 Feb 2021 09:52:15 -0500 |
parents | 6d35ab1cd544 |
children | 7854d5752dd2 |
line wrap: on
line source
######################################################################## ## ## Copyright (C) 2008-2021 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. ## ## 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 ## <https://www.gnu.org/licenses/>. ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}) ## @deftypefnx {} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol}) ## @deftypefnx {} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}) ## @deftypefnx {} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{}) ## @deftypefnx {} {[@var{q}, @var{nfun}] =} quadv (@dots{}) ## ## Numerically evaluate the integral of @var{f} from @var{a} to @var{b} ## using an adaptive Simpson's rule. ## ## @var{f} is a function handle, inline function, or string containing the name ## of the function to evaluate. @code{quadv} is a vectorized version of ## @code{quad} and the function defined by @var{f} must accept a scalar or ## vector as input and return a scalar, vector, or array as output. ## ## @var{a} and @var{b} are the lower and upper limits of integration. Both ## limits must be finite. ## ## The optional argument @var{tol} defines the absolute tolerance used to stop ## the adaptation procedure. The default value is 1e-6. ## ## The algorithm used by @code{quadv} involves recursively subdividing the ## integration interval and applying Simpson's rule on each subinterval. ## If @var{trace} is true then after computing each of these partial ## integrals display: (1) the total number of function evaluations, ## (2) the left end of the subinterval, (3) the length of the subinterval, ## (4) the approximation of the integral over the subinterval. ## ## Additional arguments @var{p1}, etc., are passed directly to the function ## @var{f}. To use default values for @var{tol} and @var{trace}, one may pass ## empty matrices ([]). ## ## The result of the integration is returned in @var{q}. ## ## The optional output @var{nfun} indicates the total number of function ## evaluations performed. ## ## Note: @code{quadv} is written in Octave's scripting language and can be ## used recursively in @code{dblquad} and @code{triplequad}, unlike the ## @code{quad} function. ## @seealso{quad, quadl, quadgk, quadcc, trapz, dblquad, triplequad, integral, ## integral2, integral3} ## @end deftypefn ## Algorithm: See https://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method ## for basic explanation. See NOTEs and FIXME for Octave modifications. function [q, nfun] = quadv (f, a, b, tol = [], trace = false, varargin) if (nargin < 3) print_usage (); endif if (isa (a, "single") || isa (b, "single")) eps = eps ("single"); else eps = eps ("double"); endif if (isempty (tol)) tol = 1e-6; elseif (! isscalar (tol) || tol < 0) error ("quadv: TOL must be a scalar >=0"); endif if (trace) ## Print column headers once above trace display. printf (" nfun a (b - a) q_interval\n"); endif ## NOTE: Split the interval into 3 parts which avoids problems with periodic ## functions when a, b, and (a + b)/2 fall on boundaries such as 0, 2*pi. ## For compatibility with Matlab, split in to two equal size regions on the ## left and right, and one larger central region. alpha = 0.27158; # factor for region 1 & region 3 size (~27%) b1 = a + alpha * (b - a); b2 = b - alpha * (b - a); c1 = (a + b1) / 2; c2 = (a + b) / 2; c3 = (b2 + b) / 2; fa = feval (f, a, varargin{:}); fc1 = feval (f, c1, varargin{:}); fb1 = feval (f, b1, varargin{:}); fc2 = feval (f, c2, varargin{:}); fb2 = feval (f, b2, varargin{:}); fc3 = feval (f, c3, varargin{:}); fb = feval (f, b, varargin{:}); nfun = 7; ## NOTE: If there are edge singularities, move edge point by eps*(b-a) as ## discussed in Shampine paper used to implement quadgk. if (any (isinf (fa(:)))) fa = feval (f, a + eps * (b-a), varargin{:}); nfun++; endif if (any (isinf (fb(:)))) fb = feval (f, b - eps * (b-a), varargin{:}); nfun++; endif ## Region 1 h = (b1 - a); q1 = h / 6 * (fa + 4*fc1 + fb1); [q1, nfun, hmin1] = simpsonstp (f, a, b1, c1, fa, fb1, fc1, q1, tol, nfun, abs (h), trace, varargin{:}); ## Region 2 h = (b2 - b1); q2 = h / 6 * (fb1 + 4*fc2 + fb2); [q2, nfun, hmin2] = simpsonstp (f, b1, b2, c2, fb1, fb2, fc2, q2, tol, nfun, abs (h), trace, varargin{:}); ## Region 3 h = (b - b2); q3 = h / 6 * (fb2 + 4*fc3 + fb); [q3, nfun, hmin3] = simpsonstp (f, b2, b, c3, fb2, fb, fc3, q3, tol, nfun, abs (h), trace, varargin{:}); ## Total integral over all 3 regions and verify results q = q1 + q2 + q3; hmin = min ([hmin1, hmin2, hmin3]); if (nfun > 10_000) warning ("quadv: maximum iteration count reached -- possible singular integral"); elseif (any (! isfinite (q(:)))) warning ("quadv: infinite or NaN function evaluations were returned"); elseif (hmin < (b - a) * eps) warning ("quadv: minimum step size reached -- possible singular integral"); endif endfunction function [q, nfun, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, q0, tol, nfun, hmin, trace, varargin) if (nfun > 10_000) # stop endless recursion q = q0; return; endif d = (a + c) / 2; e = (c + b) / 2; fd = feval (f, d, varargin{:}); fe = feval (f, e, varargin{:}); nfun += 2; q1 = (c - a) / 6 * (fa + 4*fd + fc); q2 = (b - c) / 6 * (fc + 4*fe + fb); q = q1 + q2; if (abs (a - c) < hmin) hmin = abs (a - c); endif delta = q - q0; # error term between new estimate and old estimate if (trace) printf ("%5d %#14.11g %16.10e %-16.11g\n", nfun, a, b-a, q + delta/15); endif ## NOTE: Not vectorizing q-q0 in the norm provides a more rigid criterion ## for matrix-valued functions. if (norm (delta, Inf) > 15*tol) ## FIXME: To keep sum of sub-interval integrands within overall tolerance ## each bisection interval should use tol/2. However, Matlab does not ## do this, and it would also profoundly increase the number of function ## evaluations required. [q1, nfun, hmin] = simpsonstp (f, a, c, d, fa, fc, fd, q1, tol, nfun, hmin, trace, varargin{:}); [q2, nfun, hmin] = simpsonstp (f, c, b, e, fc, fb, fe, q2, tol, nfun, hmin, trace, varargin{:}); q = q1 + q2; else q += delta / 15; # NOTE: Richardson extrapolation correction endif endfunction %!assert (quadv (@sin, 0, 2*pi), 0, 1e-6) %!assert (quadv (@sin, 0, pi), 2, 1e-6) ## Test weak singularities at the edge %!assert (quadv (@(x) 1 ./ sqrt (x), 0, 1), 2, 15*1e-6); ## Test vector-valued functions %!assert (quadv (@(x) [(sin (x)), (sin (2 * x))], 0, pi), [2, 0], 1e-6) ## Test matrix-valued functions %!assert (quadv (@(x) [ x,x,x; x,1./sqrt(x),x; x,x,x ], 0, 1), %! [0.5,0.5,0.5; 0.5,2,0.5; 0.5,0.5,0.5], 15*1e-6); ## Test periodic function %!assert <*57603> (quadv (@(t) sin (t) .^ 2, 0, 8*pi), 4*pi, 1e-6) ## Test input validation %!error quadv () %!error quadv (@sin) %!error quadv (@sin,1) %!error <TOL must be a scalar> quadv (@sin,0,1, ones (2,2)) %!error <TOL must be .* .=0> quadv (@sin,0,1, -1)