octave_packages/m/general/quadgk.m

   1 ## Copyright (C) 2008-2012 David Bateman
   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} =} quadgk (@var{f}, @var{a}, @var{b})
  21 ## @deftypefnx {Function File} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}, @var{abstol})
  22 ## @deftypefnx {Function File} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}, @var{abstol}, @var{trace})
  23 ## @deftypefnx {Function File} {@var{q} =} quadgk (@var{f}, @var{a}, @var{b}, @var{prop}, @var{val}, @dots{})
  24 ## @deftypefnx {Function File} {[@var{q}, @var{err}] =} quadgk (@dots{})
  25 ##
  26 ## Numerically evaluate the integral of @var{f} from @var{a} to @var{b}
  27 ## using adaptive Gauss-Konrod quadrature.
  28 ## @var{f} is a function handle, inline function, or string
  29 ## containing the name of the function to evaluate.
  30 ## The formulation is based on a proposal by L.F. Shampine,
  31 ## @cite{"Vectorized adaptive quadrature in @sc{matlab}", Journal of
  32 ## Computational and Applied Mathematics, pp131-140, Vol 211, Issue 2,
  33 ## Feb 2008} where all function evaluations at an iteration are
  34 ## calculated with a single call to @var{f}.  Therefore, the function
  35 ## @var{f} must be vectorized and must accept a vector of input values @var{x}
  36 ## and return an output vector representing the function evaluations at the
  37 ## given values of @var{x}.
  38 ##
  39 ## @var{a} and @var{b} are the lower and upper limits of integration.  Either
  40 ## or both limits may be infinite or contain weak end singularities.
  41 ## Variable transformation will be used to treat any infinite intervals and
  42 ## weaken the singularities.  For example:
  43 ##
  44 ## @example
  45 ## quadgk (@@(x) 1 ./ (sqrt (x) .* (x + 1)), 0, Inf)
  46 ## @end example
  47 ##
  48 ## @noindent
  49 ## Note that the formulation of the integrand uses the
  50 ## element-by-element operator @code{./} and all user functions to
  51 ## @code{quadgk} should do the same.
  52 ##
  53 ## The optional argument @var{tol} defines the absolute tolerance used to stop
  54 ## the integration procedure.  The default value is @math{1e^{-10}}.
  55 ##
  56 ## The algorithm used by @code{quadgk} involves subdividing the
  57 ## integration interval and evaluating each subinterval.
  58 ## If @var{trace} is true then after computing each of these partial
  59 ## integrals display: (1) the number of subintervals at this step,
  60 ## (2) the current estimate of the error @var{err}, (3) the current estimate
  61 ## for the integral @var{q}.
  62 ##
  63 ## Alternatively, properties of @code{quadgk} can be passed to the function as
  64 ## pairs @code{"@var{prop}", @var{val}}.  Valid properties are
  65 ##
  66 ## @table @code
  67 ## @item AbsTol
  68 ## Define the absolute error tolerance for the quadrature.  The default
  69 ## absolute tolerance is 1e-10.
  70 ##
  71 ## @item RelTol
  72 ## Define the relative error tolerance for the quadrature.  The default
  73 ## relative tolerance is 1e-5.
  74 ##
  75 ## @item MaxIntervalCount
  76 ## @code{quadgk} initially subdivides the interval on which to perform
  77 ## the quadrature into 10 intervals.  Subintervals that have an
  78 ## unacceptable error are subdivided and re-evaluated.  If the number of
  79 ## subintervals exceeds 650 subintervals at any point then a poor
  80 ## convergence is signaled and the current estimate of the integral is
  81 ## returned.  The property 'MaxIntervalCount' can be used to alter the
  82 ## number of subintervals that can exist before exiting.
  83 ##
  84 ## @item WayPoints
  85 ## Discontinuities in the first derivative of the function to integrate can be
  86 ## flagged with the  @code{"WayPoints"} property.  This forces the ends of
  87 ## a subinterval to fall on the breakpoints of the function and can result in
  88 ## significantly improved estimation of the error in the integral, faster
  89 ## computation, or both.  For example,
  90 ##
  91 ## @example
  92 ## quadgk (@@(x) abs (1 - x.^2), 0, 2, "Waypoints", 1)
  93 ## @end example
  94 ##
  95 ## @noindent
  96 ## signals the breakpoint in the integrand at @code{@var{x} = 1}.
  97 ##
  98 ## @item Trace
  99 ## If logically true @code{quadgk} prints information on the
 100 ## convergence of the quadrature at each iteration.
 101 ## @end table
 102 ##
 103 ## If any of @var{a}, @var{b}, or @var{waypoints} is complex then the
 104 ## quadrature is treated as a contour integral along a piecewise
 105 ## continuous path defined by the above.  In this case the integral is
 106 ## assumed to have no edge singularities.  For example,
 107 ##
 108 ## @example
 109 ## @group
 110 ## quadgk (@@(z) log (z), 1+1i, 1+1i, "WayPoints",
 111 ##         [1-1i, -1,-1i, -1+1i])
 112 ## @end group
 113 ## @end example
 114 ##
 115 ## @noindent
 116 ## integrates @code{log (z)} along the square defined by @code{[1+1i,
 117 ##  1-1i, -1-1i, -1+1i]}
 118 ##
 119 ## The result of the integration is returned in @var{q}.
 120 ## @var{err} is an approximate bound on the error in the integral
 121 ## @code{abs (@var{q} - @var{I})}, where @var{I} is the exact value of the
 122 ## integral.
 123 ##
 124 ## @seealso{quad, quadv, quadl, quadcc, trapz, dblquad, triplequad}
 125 ## @end deftypefn
 126
 127 function [q, err] = quadgk (f, a, b, varargin)
 128   if (nargin < 3)
 129     print_usage ();
 130   endif
 131
 132   if (b < a)
 133     [q, err] = quadgk (f, b, a, varargin{:});
 134     q = -q;
 135   else
 136     abstol = 1e-10;
 137     reltol = 1e-5;
 138     waypoints = [];
 139     maxint = 650;
 140     trace = false;
 141
 142     if (nargin > 3)
 143       if (! ischar (varargin{1}))
 144         if (!isempty (varargin{1}))
 145           abstol = varargin{1};
 146           reltol = 0;
 147         endif
 148         if (nargin > 4)
 149           trace = varargin{2};
 150         endif
 151         if (nargin > 5)
 152           error ("quadgk: can not pass additional arguments to user function");
 153         endif
 154       else
 155         idx = 1;
 156         while (idx < nargin - 3)
 157           if (ischar (varargin{idx}))
 158             str = varargin{idx++};
 159             if (strcmpi (str, "reltol"))
 160               reltol = varargin{idx++};
 161             elseif (strcmpi (str, "abstol"))
 162               abstol = varargin{idx++};
 163             elseif (strcmpi (str, "waypoints"))
 164               waypoints = varargin{idx++} (:);
 165               if (isreal(waypoints))
 166                 waypoints (waypoints < a | waypoints > b) = [];
 167               endif
 168             elseif (strcmpi (str, "maxintervalcount"))
 169               maxint = varargin{idx++};
 170             elseif (strcmpi (str, "trace"))
 171               trace = varargin{idx++};
 172             else
 173               error ("quadgk: unknown property %s", str);
 174             endif
 175           else
 176             error ("quadgk: expecting property to be a string");
 177           endif
 178         endwhile
 179         if (idx != nargin - 2)
 180           error ("quadgk: expecting properties in pairs");
 181         endif
 182       endif
 183     endif
 184
 185     ## Convert function given as a string to a function handle
 186     if (ischar (f))
 187       f = @(x) feval (f, x);
 188     endif
 189
 190     ## Use variable subsitution to weaken endpoint singularities and to
 191     ## perform integration with endpoints at infinity. No transform for
 192     ## contour integrals
 193     if (iscomplex (a) || iscomplex (b) || iscomplex(waypoints))
 194       ## contour integral, no transform
 195       subs = [a; waypoints; b];
 196       h = sum (abs (diff (subs)));
 197       h0 = h;
 198       trans = @(t) t;
 199     elseif (isinf (a) && isinf(b))
 200       ## Standard Infinite to finite integral transformation.
 201       ##   \int_{-\infinity_^\infinity f(x) dx = \int_-1^1 f (g(t)) g'(t) dt
 202       ## where
 203       ##   g(t)  = t / (1 - t^2)
 204       ##   g'(t) =  (1 + t^2) / (1 - t^2) ^ 2
 205       ## waypoint transform is then
 206       ##   t =  (2 * g(t)) ./ (1 + sqrt(1 + 4 * g(t) .^ 2))
 207       if (!isempty (waypoints))
 208         trans = @(x) (2 * x) ./ (1 + sqrt(1 + 4 * x .^ 2));
 209         subs = [-1; trans(waypoints); 1];
 210       else
 211         subs = linspace (-1, 1, 11)';
 212       endif
 213       h = 2;
 214       h0 = b - a;
 215       trans = @(t) t ./ (1 - t.^2);
 216       f = @(t) f (t ./ (1 - t .^ 2)) .* (1 + t .^ 2) ./ ((1 - t .^ 2) .^ 2);
 217     elseif (isinf(a))
 218       ## Formula defined in Shampine paper as two separate steps. One to
 219       ## weaken singularity at finite end, then a second to transform to
 220       ## a finite interval. The singularity weakening transform is
 221       ##   \int_{-\infinity}^b f(x) dx =
 222       ##               - \int_{-\infinity}^0 f (b - t^2) 2 t dt
 223       ## (note minus sign) and the finite interval transform is
 224       ##   \int_{-\infinity}^0 f(b - t^2)  2 t dt =
 225       ##                  \int_{-1}^0 f (b - g(s) ^ 2) 2 g(s) g'(s) ds
 226       ## where
 227       ##   g(s)  = s / (1 + s)
 228       ##   g'(s) = 1 / (1 + s) ^ 2
 229       ## waypoint transform is then
 230       ##   t = sqrt (b - x)
 231       ##   s =  - t / (t + 1)
 232       if (!isempty (waypoints))
 233         tmp = sqrt (b - waypoints);
 234         trans = @(x)  - x ./ (x + 1);
 235         subs = [-1; trans(tmp); 0];
 236       else
 237         subs = linspace (-1, 0, 11)';
 238       endif
 239       h = 1;
 240       h0 = b - a;
 241       trans = @(t) b - (t ./ (1 + t)).^2;
 242       f = @(s) - 2 * s .* f (b -  (s ./ (1 + s)) .^ 2) ./ ((1 + s) .^ 3);
 243     elseif (isinf(b))
 244       ## Formula defined in Shampine paper as two separate steps. One to
 245       ## weaken singularity at finite end, then a second to transform to
 246       ## a finite interval. The singularity weakening transform is
 247       ##   \int_a^\infinity f(x) dx = \int_0^\infinity f (a + t^2) 2 t dt
 248       ## and the finite interval transform is
 249       ##  \int_0^\infinity f(a + t^2)  2 t dt =
 250       ##           \int_0^1 f (a + g(s) ^ 2) 2 g(s) g'(s) ds
 251       ## where
 252       ##   g(s)  = s / (1 - s)
 253       ##   g'(s) = 1 / (1 - s) ^ 2
 254       ## waypoint transform is then
 255       ##   t = sqrt (x - a)
 256       ##   s = t / (t + 1)
 257       if (!isempty (waypoints))
 258         tmp = sqrt (waypoints - a);
 259         trans = @(x) x ./ (x + 1);
 260         subs = [0; trans(tmp); 1];
 261       else
 262         subs = linspace (0, 1, 11)';
 263       endif
 264       h = 1;
 265       h0 = b - a;
 266       trans = @(t) a + (t ./ (1 - t)).^2;
 267       f = @(s) 2 * s .* f (a +  (s ./ (1 - s)) .^ 2) ./ ((1 - s) .^ 3);
 268     else
 269       ## Davis, Rabinowitz, "Methods of Numerical Integration" p441 2ed.
 270       ## Presented in section 5 of the Shampine paper as
 271       ##   g(t) = ((b - a) / 2) * (t / 2 * (3 - t^2)) + (b + a) / 2
 272       ##   g'(t) = ((b-a)/4) * (3 - 3t^2);
 273       ## waypoint transform can then be found by solving for t with
 274       ## Maxima (solve (c + 3*t -  3^3, t);). This gives 3 roots, two of
 275       ## which are complex for values between a and b and so can be
 276       ## ignored. The third is
 277       ##  c = (-4*x + 2*(b+a)) / (b-a);
 278       ##  k = ((sqrt(c^2 - 4) + c)/2)^(1/3);
 279       ##  t = (sqrt(3)* 1i * (1 - k^2) - (1 + k^2)) / 2 / k;
 280       if (! isempty (waypoints))
 281         trans = @__quadgk_finite_waypoint__;
 282         subs = [-1; trans(waypoints, a, b); 1];
 283       else
 284         subs = linspace(-1, 1, 11)';
 285       endif
 286       h = 2;
 287       h0 = b - a;
 288       trans = @(t) ((b - a) ./ 4) * t .* (3 - t.^2) + (b + a) ./ 2;
 289       f = @(t) f((b - a) ./ 4 .* t .* (3 - t.^2) + (b + a) ./ 2) .* ...
 290            3 .* (b - a) ./ 4 .* (1 - t.^2);
 291     endif
 292
 293     ## Split interval into at least 10 subinterval with a 15 point
 294     ## Gauss-Kronrod rule giving a minimum of 150 function evaluations
 295     while (length (subs) < 11)
 296       subs = [subs' ; subs(1:end-1)' + diff(subs') ./ 2, NaN](:)(1 : end - 1);
 297     endwhile
 298     subs = [subs(1:end-1), subs(2:end)];
 299
 300     warn_state = warning ("query", "Octave:divide-by-zero");
 301
 302     unwind_protect
 303       ## Singularity will cause divide by zero warnings
 304       warning ("off", "Octave:divide-by-zero");
 305
 306       ## Initial evaluation of the integrand on the subintervals
 307       [q_subs, q_errs] = __quadgk_eval__ (f, subs);
 308       q0 = sum (q_subs);
 309       err0 = sum (q_errs);
 310
 311       if (isa (a, "single") || isa (b, "single") || isa (waypoints, "single"))
 312         myeps = eps ("single");
 313       else
 314         myeps = eps;
 315       endif
 316
 317       first = true;
 318       while (true)
 319         ## Check for subintervals that are too small. Test must be
 320         ## performed in untransformed subintervals. What is a good
 321         ## value for this test. Shampine suggests 100*eps
 322         if (any (abs (diff (trans (subs), [], 2) / h0) < 100 * myeps))
 323           q = q0;
 324           err = err0;
 325           break;
 326         endif
 327
 328         ## Quit if any evaluations are not finite (Inf or NaN)
 329         if (any (! isfinite (q_subs)))
 330           warning ("quadgk: non finite integrand encountered");
 331           q = q0;
 332           err = err0;
 333           break;
 334         endif
 335
 336         tol = max (abstol, reltol .* abs (q0));
 337
 338         ## If the global error estimate is meet exit
 339         if (err0 < tol)
 340           q = q0;
 341           err = err0;
 342           break;
 343         endif
 344
 345         ## Accept the subintervals that meet the convergence criteria
 346         idx = find (abs (q_errs) < tol .* abs(diff (subs, [], 2)) ./ h);
 347         if (first)
 348           q = sum (q_subs (idx));
 349           err = sum (q_errs(idx));
 350           first = false;
 351         else
 352           q0 = q + sum (q_subs);
 353           err0 = err + sum (q_errs);
 354           q += sum (q_subs (idx));
 355           err += sum (q_errs(idx));
 356         endif
 357         subs(idx,:) = [];
 358
 359         ## If no remaining subintervals exit
 360         if (rows (subs) == 0)
 361           break;
 362         endif
 363
 364         if (trace)
 365           disp([rows(subs), err, q0]);
 366         endif
 367
 368         ## Split remaining subintervals in two
 369         mid = (subs(:,2) + subs(:,1)) ./ 2;
 370         subs = [subs(:,1), mid; mid, subs(:,2)];
 371
 372         ## If the maximum subinterval count is met accept remaining
 373         ## subinterval and exit
 374         if (rows (subs) > maxint)
 375           warning ("quadgk: maximum interval count (%d) met", maxint);
 376           q += sum (q_subs);
 377           err += sum (q_errs);
 378           break;
 379         endif
 380
 381         ## Evaluation of the integrand on the remaining subintervals
 382         [q_subs, q_errs] = __quadgk_eval__ (f, subs);
 383       endwhile
 384
 385       if (err > max (abstol, reltol * abs(q)))
 386         warning ("quadgk: Error tolerance not met. Estimated error %g", err);
 387       endif
 388     unwind_protect_cleanup
 389       if (strcmp (warn_state.state, "on"))
 390         warning ("on", "Octave:divide-by-zero");
 391       endif
 392     end_unwind_protect
 393   endif
 394 endfunction
 395
 396 function [q, err] = __quadgk_eval__ (f, subs)
 397   ## A (15,7) point pair of Gauss-Konrod quadrature rules. The abscissa
 398   ## and weights are copied directly from dqk15w.f from quadpack
 399
 400   persistent abscissa = [-0.9914553711208126e+00, -0.9491079123427585e+00, ...
 401                          -0.8648644233597691e+00, -0.7415311855993944e+00, ...
 402                          -0.5860872354676911e+00, -0.4058451513773972e+00, ...
 403                          -0.2077849550078985e+00,  0.0000000000000000e+00, ...
 404                           0.2077849550078985e+00,  0.4058451513773972e+00, ...
 405                           0.5860872354676911e+00,  0.7415311855993944e+00, ...
 406                           0.8648644233597691e+00,  0.9491079123427585e+00, ...
 407                           0.9914553711208126e+00];
 408
 409   persistent weights15 = ...
 410       diag ([0.2293532201052922e-01,  0.6309209262997855e-01, ...
 411              0.1047900103222502e+00,  0.1406532597155259e+00, ...
 412              0.1690047266392679e+00,  0.1903505780647854e+00, ...
 413              0.2044329400752989e+00,  0.2094821410847278e+00, ...
 414              0.2044329400752989e+00,  0.1903505780647854e+00, ...
 415              0.1690047266392679e+00,  0.1406532597155259e+00, ...
 416              0.1047900103222502e+00,  0.6309209262997855e-01, ...
 417              0.2293532201052922e-01]);
 418
 419   persistent weights7  = ...
 420       diag ([0.1294849661688697e+00,  0.2797053914892767e+00, ...
 421              0.3818300505051889e+00,  0.4179591836734694e+00, ...
 422              0.3818300505051889e+00,  0.2797053914892767e+00, ...
 423              0.1294849661688697e+00]);
 424
 425   halfwidth = diff (subs, [], 2) ./ 2;
 426   center = sum (subs, 2) ./ 2;;
 427   x = bsxfun (@plus, halfwidth * abscissa, center);
 428   y = reshape (f (x(:)), size(x));
 429
 430   ## This is faster than using bsxfun as the * operator can use a
 431   ## single BLAS call, rather than rows(sub) calls to the @times
 432   ## function.
 433   q = sum (y * weights15, 2) .* halfwidth;
 434   err = abs (sum (y(:,2:2:end) * weights7, 2) .* halfwidth - q);
 435 endfunction
 436
 437 function t = __quadgk_finite_waypoint__ (x, a, b)
 438   c = (-4 .* x + 2.* (b + a)) ./ (b - a);
 439   k = ((sqrt(c .^ 2 - 4) + c) ./ 2) .^ (1/3);
 440   t = real ((sqrt(3) .* 1i * (1 - k .^ 2) - (1 + k .^ 2)) ./ 2 ./ k);
 441 endfunction
 442
 443 %error (quadgk (@sin))
 444 %error (quadgk (@sin, -pi))
 445 %error (quadgk (@sin, -pi, pi, 'DummyArg'))
 446
 447 %!assert (quadgk(@sin,-pi,pi), 0, 1e-6)
 448 %!assert (quadgk(inline('sin'),-pi,pi), 0, 1e-6)
 449 %!assert (quadgk('sin',-pi,pi), 0, 1e-6)
 450 %!assert (quadgk(@sin,-pi,pi,'waypoints', 0, 'MaxIntervalCount', 100, 'reltol', 1e-3, 'abstol', 1e-6, 'trace', false), 0, 1e-6)
 451 %!assert (quadgk(@sin,-pi,pi,1e-6,false), 0, 1e-6)
 452
 453 %!assert (quadgk(@sin,-pi,0), -2, 1e-6)
 454 %!assert (quadgk(@sin,0,pi), 2, 1e-6)
 455 %!assert (quadgk(@(x) 1./sqrt(x), 0, 1), 2, 1e-6)
 456 %!assert (quadgk (@(x) abs (1 - x.^2), 0, 2, 'Waypoints', 1), 2, 1e-6)
 457 %!assert (quadgk(@(x) 1./(sqrt(x).*(x+1)), 0, Inf), pi, 1e-6)
 458 %!assert (quadgk (@(z) log (z), 1+1i, 1+1i, 'WayPoints', [1-1i, -1,-1i, -1+1i]), -pi * 1i, 1e-6)
 459
 460 %!assert (quadgk (@(x) exp(-x .^ 2), -Inf, Inf), sqrt(pi), 1e-6)
 461 %!assert (quadgk (@(x) exp(-x .^ 2), -Inf, 0), sqrt(pi)/2, 1e-6)