octave_packages/m/general/quadl.m

   1 ## Copyright (C) 1998-2012 Walter Gautschi
   2 ##
   3 ## This file is part of Octave.
   4 ##
   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.
   9 ##
  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.
  14 ##
  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/>.
  18
  19 ## -*- texinfo -*-
  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{})
  24 ##
  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.
  31 ##
  32 ## @var{a} and @var{b} are the lower and upper limits of integration.  Both
  33 ## limits must be finite.
  34 ##
  35 ## The optional argument @var{tol} defines the relative tolerance with which
  36 ## to perform the integration.  The default value is @code{eps}.
  37 ##
  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.
  43 ##
  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 ([]).
  47 ##
  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}
  52 ## @end deftypefn
  53
  54 ##   Author: Walter Gautschi
  55 ##   Date: 08/03/98
  56 ##   Reference: Gander, Computermathematik, Birkhaeuser, 1992.
  57
  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
  64
  65 function q = quadl (f, a, b, tol = [], trace = false, varargin)
  66
  67   if (nargin < 3)
  68     print_usage ();
  69   endif
  70
  71   if (isa (a, "single") || isa (b, "single"))
  72     myeps = eps ("single");
  73   else
  74     myeps = eps;
  75   endif
  76   if (isempty (tol))
  77     tol = myeps;
  78   endif
  79   if (isempty (trace))
  80     trace = false;
  81   endif
  82   if (tol < myeps)
  83     tol = myeps;
  84   endif
  85
  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;
  92
  93   m = (a+b)/2;
  94   h = (b-a)/2;
  95   alpha = sqrt (2/3);
  96   beta = 1/sqrt (5);
  97
  98   x1 = .942882415695480;
  99   x2 = .641853342345781;
 100   x3 = .236383199662150;
 101
 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];
 104
 105   y = feval (f, x, varargin{:});
 106
 107   fa = y(1);
 108   fb = y(13);
 109
 110   i2 = (h/6)*(y(1) + y(13) + 5*(y(5)+y(9)));
 111
 112   i1 = (h/1470)*(   77*(y(1)+y(13))
 113                  + 432*(y(3)+y(11))
 114                  + 625*(y(5)+y(9))
 115                  + 672*y(7));
 116
 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));
 124
 125   s = sign (is);
 126   if (s == 0)
 127     s = 1;
 128   endif
 129   erri1 = abs (i1-is);
 130   erri2 = abs (i2-is);
 131   if (erri2 != 0)
 132     R = erri1/erri2;
 133   else
 134     R = 1;
 135   endif
 136   if (R > 0 && R < 1)
 137     tol = tol/R;
 138   endif
 139   is = s * abs(is) * tol/myeps;
 140   if (is == 0)
 141     is = b-a;
 142   endif
 143
 144   q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin{:});
 145
 146 endfunction
 147
 148 ## ADAPTLOBSTP  Recursive function used by QUADL.
 149 ##
 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.
 155 ##
 156 ##   Walter Gautschi, 08/03/98
 157
 158 function q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin)
 159   global __quadl_recurse_done__;
 160   global __quadl_need_warning__;
 161
 162   h = (b-a)/2;
 163   m = (a+b)/2;
 164   alpha = sqrt (2/3);
 165   beta = 1 / sqrt(5);
 166   mll = m-alpha*h;
 167   ml  = m-beta*h;
 168   mr  = m+beta*h;
 169   mrr = m+alpha*h;
 170   x = [mll, ml, m, mr, mrr];
 171   y = feval (f, x, varargin{:});
 172   fmll = y(1);
 173   fml  = y(2);
 174   fm   = y(3);
 175   fmr  = y(4);
 176   fmrr = y(5);
 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;
 184     endif
 185     q = i1;
 186     if (trace)
 187       disp ([a, b-a, q]);
 188     endif
 189   else
 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{:}));
 197   endif
 198 endfunction
 199
 200
 201 ## basic functionality
 202 %!assert (quadl (@(x) sin (x), 0, pi, [], []), 2, -3e-16)
 203
 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)
 207
 208 ## extra parameters
 209 %!assert (quadl (@(x,a,b) sin (a + b*x), 0, 1, [], [], 2, 3),
 210 %!        cos(2)/3 - cos(5)/3, -3e-16)
 211
 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)
 217