1 ## Copyright (C) 1998-2012 Walter Gautschi
3 ## This file is part of Octave.
5 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by
7 ## the Free Software Foundation; either version 3 of the License, or (at
8 ## your option) any later version.
10 ## Octave is distributed in the hope that it will be useful, but
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 ## General Public License for more details.
15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>.
20 ## @deftypefn {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b})
21 ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol})
22 ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace})
23 ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{})
25 ## Numerically evaluate the integral of @var{f} from @var{a} to @var{b}
26 ## using an adaptive Lobatto rule.
27 ## @var{f} is a function handle, inline function, or string
28 ## containing the name of the function to evaluate.
29 ## The function @var{f} must be vectorized and return a vector of output values
30 ## if given a vector of input values.
32 ## @var{a} and @var{b} are the lower and upper limits of integration. Both
33 ## limits must be finite.
35 ## The optional argument @var{tol} defines the relative tolerance with which
36 ## to perform the integration. The default value is @code{eps}.
38 ## The algorithm used by @code{quadl} involves recursively subdividing the
39 ## integration interval.
40 ## If @var{trace} is defined then for each subinterval display: (1) the left
41 ## end of the subinterval, (2) the length of the subinterval, (3) the
42 ## approximation of the integral over the subinterval.
44 ## Additional arguments @var{p1}, etc., are passed directly to the function
45 ## @var{f}. To use default values for @var{tol} and @var{trace}, one may pass
46 ## empty matrices ([]).
48 ## Reference: W. Gander and W. Gautschi, @cite{Adaptive Quadrature -
49 ## Revisited}, BIT Vol. 40, No. 1, March 2000, pp. 84--101.
50 ## @url{http://www.inf.ethz.ch/personal/gander/}
51 ## @seealso{quad, quadv, quadgk, quadcc, trapz, dblquad, triplequad}
54 ## Author: Walter Gautschi
56 ## Reference: Gander, Computermathematik, Birkhaeuser, 1992.
58 ## 2003-08-05 Shai Ayal
59 ## * permission from author to release as GPL
60 ## 2004-02-10 Paul Kienzle
61 ## * renamed to quadl for compatibility
62 ## * replace global variable terminate2 with local function need_warning
63 ## * add paper ref to docs
65 function q = quadl (f, a, b, tol = [], trace = false, varargin)
71 if (isa (a, "single") || isa (b, "single"))
72 myeps = eps ("single");
86 ## Track whether recursion has occurred
87 global __quadl_recurse_done__;
88 __quadl_recurse_done__ = false;
89 ## Track whether warning about machine precision has been issued
90 global __quadl_need_warning__;
91 __quadl_need_warning__ = true;
98 x1 = .942882415695480;
99 x2 = .641853342345781;
100 x3 = .236383199662150;
102 x = [a, m-x1*h, m-alpha*h, m-x2*h, m-beta*h, m-x3*h, m, m+x3*h, ...
103 m+beta*h, m+x2*h, m+alpha*h, m+x1*h, b];
105 y = feval (f, x, varargin{:});
110 i2 = (h/6)*(y(1) + y(13) + 5*(y(5)+y(9)));
112 i1 = (h/1470)*( 77*(y(1)+y(13))
117 is = h*( .0158271919734802*(y(1)+y(13))
118 +.0942738402188500*(y(2)+y(12))
119 + .155071987336585*(y(3)+y(11))
120 + .188821573960182*(y(4)+y(10))
121 + .199773405226859*(y(5)+y(9))
122 + .224926465333340*(y(6)+y(8))
123 + .242611071901408*y(7));
139 is = s * abs(is) * tol/myeps;
144 q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin{:});
148 ## ADAPTLOBSTP Recursive function used by QUADL.
150 ## Q = ADAPTLOBSTP('F', A, B, FA, FB, IS, TRACE) tries to
151 ## approximate the integral of F(X) from A to B to
152 ## an appropriate relative error. The argument 'F' is
153 ## a string containing the name of f. The remaining
154 ## arguments are generated by ADAPTLOB or by recursion.
156 ## Walter Gautschi, 08/03/98
158 function q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin)
159 global __quadl_recurse_done__;
160 global __quadl_need_warning__;
170 x = [mll, ml, m, mr, mrr];
171 y = feval (f, x, varargin{:});
177 i2 = (h/6)*(fa + fb + 5*(fml+fmr));
178 i1 = (h/1470)*(77*(fa+fb) + 432*(fmll+fmrr) + 625*(fml+fmr) + 672*fm);
179 if ((is+(i1-i2) == is || mll <= a || b <= mrr) && __quadl_recurse_done__)
180 if ((m <= a || b <= m) && __quadl_need_warning__)
181 warning ("quadl: interval contains no more machine number");
182 warning ("quadl: required tolerance may not be met");
183 __quadl_need_warning__ = false;
190 __quadl_recurse_done__ = true;
191 q = ( adaptlobstp (f, a , mll, fa , fmll, is, trace, varargin{:})
192 + adaptlobstp (f, mll, ml , fmll, fml , is, trace, varargin{:})
193 + adaptlobstp (f, ml , m , fml , fm , is, trace, varargin{:})
194 + adaptlobstp (f, m , mr , fm , fmr , is, trace, varargin{:})
195 + adaptlobstp (f, mr , mrr, fmr , fmrr, is, trace, varargin{:})
196 + adaptlobstp (f, mrr, b , fmrr, fb , is, trace, varargin{:}));
201 ## basic functionality
202 %!assert (quadl (@(x) sin (x), 0, pi, [], []), 2, -3e-16)
204 ## the values here are very high so it may be unavoidable that this fails
205 %!assert (quadl (@(x) sin (3*x).*cosh (x).*sinh (x),10,15),
206 %! 2.588424538641647e+10, -1.1e-14)
209 %!assert (quadl (@(x,a,b) sin (a + b*x), 0, 1, [], [], 2, 3),
210 %! cos(2)/3 - cos(5)/3, -3e-16)
212 ## test different tolerances.
213 %!assert (quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.3, []),
214 %! (60 + sin(4) - sin(64))/12, -0.3)
215 %!assert (quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.1, []),
216 %! (60 + sin(4) - sin(64))/12, -0.1)