octave_packages/m/general/quadv.m

   1 ## Copyright (C) 2008-2012 David Bateman
   2 ## Copyright (C) 2012 Alexander Klein
   3 ##
   4 ## This file is part of Octave.
   5 ##
   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.
  10 ##
  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.
  15 ##
  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/>.
  19
  20 ## -*- texinfo -*-
  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{})
  26 ##
  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.
  34 ##
  35 ## @var{a} and @var{b} are the lower and upper limits of integration.  Both
  36 ## limits must be finite.
  37 ##
  38 ## The optional argument @var{tol} defines the tolerance used to stop
  39 ## the adaptation procedure.  The default value is @math{1e^{-6}}.
  40 ##
  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.
  47 ##
  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 ([]).
  51 ##
  52 ## The result of the integration is returned in @var{q}.  @var{nfun} indicates
  53 ## the number of function evaluations that were made.
  54 ##
  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}
  59 ## @end deftypefn
  60
  61 function [q, nfun] = quadv (f, a, b, tol, trace, varargin)
  62   ## TODO: Make norm for convergence testing configurable
  63
  64   if (nargin < 3)
  65     print_usage ();
  66   endif
  67   if (nargin < 4)
  68     tol = [];
  69   endif
  70   if (nargin < 5)
  71     trace = [];
  72   endif
  73   if (isa (a, "single") || isa (b, "single"))
  74     myeps = eps ("single");
  75   else
  76     myeps = eps;
  77   endif
  78   if (isempty (tol))
  79     tol = 1e-6;
  80   endif
  81   if (isempty (trace))
  82     trace = 0;
  83   endif
  84
  85   ## Split the interval into 3 abscissa, and apply a 3 point Simpson's rule
  86   c = (a + b) / 2;
  87   fa = feval (f, a, varargin{:});
  88   fc = feval (f, c, varargin{:});
  89   fb = feval (f, b, varargin{:});
  90   nfun = 3;
  91
  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{:});
  96   endif
  97   if (any (isinf (fb(:))))
  98     fb = feval (f, b - myeps * (b-a), varargin{:});
  99   endif
 100
 101   h = (b - a);
 102   q = (b - a) / 6 * (fa + 4 * fc + fb);
 103
 104   [q, nfun, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, q, nfun, abs (h),
 105                                 tol, trace, varargin{:});
 106
 107   if (nfun > 10000)
 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");
 113   endif
 114 endfunction
 115
 116 function [q, nfun, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, q0,
 117                                        nfun, hmin, tol, trace, varargin)
 118   if (nfun > 10000)
 119     q = q0;
 120   else
 121     d = (a + c) / 2;
 122     e = (c + b) / 2;
 123     fd = feval (f, d, varargin{:});
 124     fe = feval (f, e, varargin{:});
 125     nfun += 2;
 126     q1 = (c - a) / 6 * (fa + 4 * fd + fc);
 127     q2 = (b - c) / 6 * (fc + 4 * fe + fb);
 128     q = q1 + q2;
 129
 130     if (abs(a -  c) < hmin)
 131       hmin = abs (a - c);
 132     endif
 133
 134     if (trace)
 135       disp ([nfun, a, b-a, q]);
 136     endif
 137
 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{:});
 146       q = q1 + q2;
 147     endif
 148   endif
 149 endfunction
 150
 151 %!assert (quadv (@sin, 0, 2 * pi), 0, 1e-5)
 152 %!assert (quadv (@sin, 0, pi), 2, 1e-5)
 153
 154 %% Handles weak singularities at the edge
 155 %!assert (quadv (@(x) 1 ./ sqrt(x), 0, 1), 2, 1e-5)
 156
 157 %% Handles vector-valued functions
 158 %!assert (quadv (@(x) [(sin (x)), (sin (2 * x))], 0, pi), [2, 0], 1e-5)
 159
 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)