1 ## Copyright (C) 2008-2012 David Bateman
2 ## Copyright (C) 2012 Alexander Klein
4 ## This file is part of Octave.
6 ## Octave is free software; you can redistribute it and/or modify it
7 ## under the terms of the GNU General Public License as published by
8 ## the Free Software Foundation; either version 3 of the License, or (at
9 ## your option) any later version.
11 ## Octave is distributed in the hope that it will be useful, but
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 ## General Public License for more details.
16 ## You should have received a copy of the GNU General Public License
17 ## along with Octave; see the file COPYING. If not, see
18 ## <http://www.gnu.org/licenses/>.
21 ## @deftypefn {Function File} {@var{q} =} quadv (@var{f}, @var{a}, @var{b})
22 ## @deftypefnx {Function File} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol})
23 ## @deftypefnx {Function File} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace})
24 ## @deftypefnx {Function File} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{})
25 ## @deftypefnx {Function File} {[@var{q}, @var{nfun}] =} quadv (@dots{})
27 ## Numerically evaluate the integral of @var{f} from @var{a} to @var{b}
28 ## using an adaptive Simpson's rule.
29 ## @var{f} is a function handle, inline function, or string
30 ## containing the name of the function to evaluate.
31 ## @code{quadv} is a vectorized version of @code{quad} and the function
32 ## defined by @var{f} must accept a scalar or vector as input and return a
33 ## scalar, vector, or array as output.
35 ## @var{a} and @var{b} are the lower and upper limits of integration. Both
36 ## limits must be finite.
38 ## The optional argument @var{tol} defines the tolerance used to stop
39 ## the adaptation procedure. The default value is @math{1e^{-6}}.
41 ## The algorithm used by @code{quadv} involves recursively subdividing the
42 ## integration interval and applying Simpson's rule on each subinterval.
43 ## If @var{trace} is true then after computing each of these partial
44 ## integrals display: (1) the total number of function evaluations,
45 ## (2) the left end of the subinterval, (3) the length of the subinterval,
46 ## (4) the approximation of the integral over the subinterval.
48 ## Additional arguments @var{p1}, etc., are passed directly to the function
49 ## @var{f}. To use default values for @var{tol} and @var{trace}, one may pass
50 ## empty matrices ([]).
52 ## The result of the integration is returned in @var{q}. @var{nfun} indicates
53 ## the number of function evaluations that were made.
55 ## Note: @code{quadv} is written in Octave's scripting language and can be
56 ## used recursively in @code{dblquad} and @code{triplequad}, unlike the
57 ## similar @code{quad} function.
58 ## @seealso{quad, quadl, quadgk, quadcc, trapz, dblquad, triplequad}
61 function [q, nfun] = quadv (f, a, b, tol, trace, varargin)
62 ## TODO: Make norm for convergence testing configurable
73 if (isa (a, "single") || isa (b, "single"))
74 myeps = eps ("single");
85 ## Split the interval into 3 abscissa, and apply a 3 point Simpson's rule
87 fa = feval (f, a, varargin{:});
88 fc = feval (f, c, varargin{:});
89 fb = feval (f, b, varargin{:});
92 ## If have edge singularities, move edge point by eps*(b-a) as
93 ## discussed in Shampine paper used to implement quadgk
94 if (any (isinf (fa(:))))
95 fa = feval (f, a + myeps * (b-a), varargin{:});
97 if (any (isinf (fb(:))))
98 fb = feval (f, b - myeps * (b-a), varargin{:});
102 q = (b - a) / 6 * (fa + 4 * fc + fb);
104 [q, nfun, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, q, nfun, abs (h),
105 tol, trace, varargin{:});
108 warning ("maximum iteration count reached");
109 elseif (any (isnan (q)(:) | isinf (q)(:)))
110 warning ("infinite or NaN function evaluations were returned");
111 elseif (hmin < (b - a) * myeps)
112 warning ("minimum step size reached -- possibly singular integral");
116 function [q, nfun, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, q0,
117 nfun, hmin, tol, trace, varargin)
123 fd = feval (f, d, varargin{:});
124 fe = feval (f, e, varargin{:});
126 q1 = (c - a) / 6 * (fa + 4 * fd + fc);
127 q2 = (b - c) / 6 * (fc + 4 * fe + fb);
130 if (abs(a - c) < hmin)
135 disp ([nfun, a, b-a, q]);
138 ## Force at least one adpative step.
139 ## Not vectorizing q-q0 in the norm provides a more rigid criterion for
140 ## matrix-valued functions.
141 if (nfun == 5 || norm (q - q0, Inf) > tol)
142 [q1, nfun, hmin] = simpsonstp (f, a, c, d, fa, fc, fd, q1, nfun, hmin,
143 tol, trace, varargin{:});
144 [q2, nfun, hmin] = simpsonstp (f, c, b, e, fc, fb, fe, q2, nfun, hmin,
145 tol, trace, varargin{:});
151 %!assert (quadv (@sin, 0, 2 * pi), 0, 1e-5)
152 %!assert (quadv (@sin, 0, pi), 2, 1e-5)
154 %% Handles weak singularities at the edge
155 %!assert (quadv (@(x) 1 ./ sqrt(x), 0, 1), 2, 1e-5)
157 %% Handles vector-valued functions
158 %!assert (quadv (@(x) [(sin (x)), (sin (2 * x))], 0, pi), [2, 0], 1e-5)
160 %% Handles matrix-valued functions
161 %!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], 1e-5)