[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Changeset]: Add additional quadrature functions
|
From: |
David Bateman |
|
Subject: |
[Changeset]: Add additional quadrature functions |
|
Date: |
Sat, 10 May 2008 22:27:31 +0200 |
|
User-agent: |
Thunderbird 2.0.0.12 (X11/20080306) |
The paper
L.F. Shampine, "Vectorized adaptive quadrature in MATLAB", Journal of
Computational and Applied Mathematics, pp131-140, Vol 211, Issue 2, Feb 2008
describes what is contained in the Matlab quadgk function. It also gives
some details on the implementation of quadv in particular quadv's
handling of edge singularities. The information here is sufficient to
implement both of these functions.
The advantage of quadgk is that at each iteration it has a single call
to the function defining the integrand with a vector of the abscissa.
This is as opposed to quad or quadl that has multiple single abscissa
call to the integrand. Therefore, in most cases quadgk is in fact faster
than quad, even though its implemented entirely as an m-file, though
some further optimizations might be possible.
Also quadgk handles edge singularities and infinity intervals through
the use of variable substitution. So
quadgk (@(x) 1 ./ sqrt (x), 0, 1)
works fine where it'll fail with quadl currently. Adding dblquad and
triplequad is also relatively simple and included in this patch. I used
quadgk as the default quadrature for these two functions.
Changset attached
D.
# HG changeset patch
# User David Bateman <address@hidden>
# Date 1210449340 -7200
# Node ID 37ac5052ef06e755f09ccff22ed7a48f00261b99
# Parent 5f4fede2376e15e60fd5d7686027122cd204f41f
Add quadv, quadgk, dblquad and triplequad functions
diff --git a/scripts/general/Makefile.in b/scripts/general/Makefile.in
--- a/scripts/general/Makefile.in
+++ b/scripts/general/Makefile.in
@@ -36,16 +36,16 @@ SOURCES = __isequal__.m __splinen__.m ac
SOURCES = __isequal__.m __splinen__.m accumarray.m arrayfun.m \
bicubic.m bitcmp.m bitget.m bitset.m blkdiag.m cart2pol.m \
cart2sph.m cell2mat.m celldisp.m circshift.m common_size.m \
- cplxpair.m cumtrapz.m deal.m del2.m diff.m flipdim.m fliplr.m \
+ cplxpair.m cumtrapz.m dblquad.m deal.m del2.m diff.m flipdim.m fliplr.m \
flipud.m genvarname.m gradient.m idivide.m ind2sub.m int2str.m \
interp1.m interp2.m interp3.m interpn.m interpft.m \
is_duplicate_entry.m isa.m isdefinite.m isdir.m isequal.m \
isequalwithequalnans.m isscalar.m issquare.m issymmetric.m \
isvector.m logical.m logspace.m mod.m nargchk.m \
nargoutchk.m nextpow2.m nthroot.m num2str.m perror.m pol2cart.m \
- polyarea.m postpad.m prepad.m quadl.m randperm.m rat.m rem.m \
- repmat.m rot90.m rotdim.m runlength.m shift.m shiftdim.m sortrows.m \
- sph2cart.m strerror.m structfun.m sub2ind.m trapz.m tril.m triu.m
+ polyarea.m postpad.m prepad.m quadgk.m quadl.m quadv.m randperm.m rat.m \
+ rem.m repmat.m rot90.m rotdim.m runlength.m shift.m shiftdim.m sortrows.m \
+ sph2cart.m strerror.m structfun.m sub2ind.m triplequad.m trapz.m tril.m
triu.m
DISTFILES = $(addprefix $(srcdir)/, Makefile.in $(SOURCES))
diff --git a/scripts/general/dblquad.m b/scripts/general/dblquad.m
new file mode 100644
--- /dev/null
+++ b/scripts/general/dblquad.m
@@ -0,0 +1,68 @@
+## Copyright (C) 2008 David Bateman
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} dblquad (@var{f}, @var{xa}, @var{xb},
@var{ya}, @var{yb}, @var{tol}, @var{quadf}, @dots{})
+## Numerically evaluate a double integral. The function over with to
+## integrate is defined by @address@hidden, and the interval for the
+## integration is defined by @address@hidden, @var{xb}, @var{ya},
+## @var{yb}]}. The function @var{f} must accept a vector @var{x} and a
+## scalar @var{y}, and return a vector of the same length as @var{x}.
+##
+## If defined, @var{tol} defines the absolute tolerance to which to
+## which to integrate each sub-integral.
+##
+## Additional arguments, are passed directly to @var{f}. To use the default
+## value for @var{tol} one may pass an empty matrix.
+## @seealso{triplequad, quad, quadv, quadl, quadgk, trapz}
+## @end deftypefn
+
+function Q = dblquad(f, xa, xb, ya, yb, tol, quadf, varargin)
+ if (nargin < 5)
+ print_usage ();
+ endif
+ if (nargin < 6 || isempty (tol))
+ tol = 1e-6;
+ endif
+ if (nargin < 7 || isempty (quadf))
+ quadf = @quadgk;
+ endif
+
+ inner = @__dblquad_inner__;
+ if (ischar (f))
+ f = @(x,y) feval (f, x, y, varargin{:});
+ varargin = {};
+ endif
+
+ Q = feval (quadf, @(y) inner (y, f, xa, xb, tol, quadf,
+ varargin{:}), ya, yb, tol);
+endfunction
+
+function Q = __dblquad_inner__ (y, f, xa, xb, tol, quadf, varargin)
+ Q = zeros (size(y));
+ for i = 1 : length (y)
+ Q(i) = feval (quadf, @(x) f(x, y(i), varargin{:}), xa, xb, tol);
+ endfor
+endfunction
+
+%% Nasty integrand to show quadgk off
+%!assert (dblquad (@(x,y) 1 ./ (x+y), 0, 1, 0, 1), 2*log(2), 1e-6)
+
+%!assert (dblquad (@(x,y) exp(-x.^2 - y.^2) , -1, 1, -1, 1, [], @quadgk), pi
* erf(1).^2, 1e-6)
+%!assert (dblquad (@(x,y) exp(-x.^2 - y.^2) , -1, 1, -1, 1, [], @quadl), pi *
erf(1).^2, 1e-6)
+%!assert (dblquad (@(x,y) exp(-x.^2 - y.^2) , -1, 1, -1, 1, [], @quadv), pi *
erf(1).^2, 1e-6)
diff --git a/scripts/general/quadgk.m b/scripts/general/quadgk.m
new file mode 100644
--- /dev/null
+++ b/scripts/general/quadgk.m
@@ -0,0 +1,435 @@
+## Copyright (C) 2008 David Bateman
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} quadgk (@var{f}, @var{a}, @var{b},
@var{abstol}, @var{trace})
+## @deftypefnx {Function File} {} quadgk (@var{f}, @var{a}, @var{b},
@var{prop}, @var{val}, @dots{})
+## @deftypefnx {Function File} address@hidden, @var{err}] =} quadgk (@dots{})
+## Numerically evaluate integral using adaptive Guass-Konrod quadrature.
+## The formulation is based on a proposal by L.F. Shampine,
+## @cite{"Vectorized adaptive quadrature in MATLAB", Journal of
+## Computational and Applied Mathematics, pp131-140, Vol 211, Issue 2,
+## Feb 2008} where all function evalutions at an iteration are
+## calculated with a single call to @var{f}. Therefore the function
+## @var{f} must be of the form @address@hidden (@var{x})} and accept
+## vector values of @var{x} and return a vector of the same length
+## representing the function evalutaions at the given values of @var{x}.
+## The function @var{f} can be defined in terms of a function handle,
+## inline function or string.
+##
+## The bounds of the quadrature @address@hidden, @var{b}] can be finite
+## or infinite and contain weak end singularities. Variable
+## transformation will be used to treat infinite intervals and weaken
+## the singularities. For example
+##
+## @example
+## quadgk(@@(x) 1 ./ (sqrt (x) .* (x + 1)), 0, Inf)
+## @end example
+##
+## @noindent
+## Note that the formulation of the integrand uses the
+## element-by-element operator @code{./} and all user functions to
+## @code{quadgk} should do the same.
+##
+## The absolute tolerance can be passed as a fourth argument in a manner
+## compatible with @code{quadv}. Equally the user can request that
+## information on the convergence can be printed is the fifth argument
+## is logicallly true.
+##
+## Alternatively, certain properties of @code{quadgk} can be passed as
+## pairs @address@hidden, @var{val}}. Valid properties are
+##
+## @table @code
+## @item AbsTol
+## Defines the absolute error tolerance for the quadrature. The default
+## absolute tolerance is 1e-10.
+##
+## @item RelTol
+## Defines the relative error tolerance for the quadrature. The default
+## relative tolerance is 1e-5.
+##
+## @item MaxIntervalCount
+## @code{quadgk} initially subdivides the interval on which to perform
+## the quadrature into 10 intervals. Sub-intervals that have an
+## unacceptable error are sub-divided and re-evaluated. If the number of
+## sub-intervals exceeds at any point 650 sub-intervals then a poor
+## convergence is signaled and the current estimate of the integral is
+## returned. The property 'MaxIntervalCount' can be used to alter the
+## number of sub-intervals that can exist before exiting.
+##
+## @item WayPoints
+## If there exists discontinuities in the first derivative of the
+## function to integrate, then these can be flagged with teh 'WayPoints'
+## property. This forces the ends of a sub-interval to fall on the
+## breakpoints of the function and can result in significantly improved
+## estimated of the error in the integral, faster computation or both.
+## For example
+##
+## @example
+## quadgk (@@(x) abs (1 - x .^ 2), 0, 2, 'Waypoints', 1)
+## @end example
+##
+## @noindent
+## signals the breakpoint in the integrand at @address@hidden = 1}.
+##
+## @item Trace
+## If logically true, then @code{quadgk} prints information on the
+## convergence of the quadrature at each iteration.
address@hidden table
+##
+## If any of @var{a}, @var{b} or @var{waypoints} is complex, then the
+## quadrature is treated as a contour integral along a piecewise
+## continuous path defined by the above. In this case the integral is
+## assuemd to have no edge singularities. For example
+##
+## @example
+## quadgk (@@ (z) log (z), 1+1i, 1+1i, 'WayPoints', [1-1i, -1,-1i, -1+1i])
+## @end example
+##
+## @noindent
+## integrates @code{log (z)} along the square defined by @code{[1+1i,
+## 1-1i, -1-1i, -1+1i]}
+##
+## If two output arguments are requested, then @var{err} returns the
+## approximate bounds on the error in the integral @code{abs (@var{q} -
+## @var{i})}, where @var{i} is the exact value of the integral.
+##
+## @seealso{triplequad, dblquad, quad, quadl, quadv, trapz}
+## @end deftypefn
+
+function [q, err] = quadgk (f, a, b, varargin)
+ if (nargin < 3)
+ print_usage ();
+ endif
+
+ if (b < a)
+ [q, err] = quadgk (f, b, a, varargin{:});
+ q = -q;
+ else
+ abstol = 1e-10;
+ reltol = 1e-5;
+ waypoints = [];
+ maxint = 650;
+ trace = false;
+
+ if (nargin > 3)
+ if (! ischar (varargin{1}))
+ if (!isempty (varargin{1}))
+ abstol = varargin{1};
+ reltol = 0;
+ endif
+ if (nargin > 4)
+ trace = varargin{2};
+ endif
+ if (nargin > 5)
+ error ("quadgk: can not pass additional arguments to user function");
+ endif
+ else
+ idx = 1;
+ while (idx < nargin - 3)
+ if (ischar (varargin{idx}))
+ str = varargin{idx++};
+ if (strcmpi (str, "reltol"))
+ reltol = varargin{idx++};
+ elseif (strcmpi (str, "abstol"))
+ abstol = varargin{idx++};
+ elseif (strcmpi (str, "waypoints"))
+ waypoints = varargin{idx++} (:);
+ if (isreal(waypoints))
+ waypoints (waypoints < a | waypoints > b) = [];
+ endif
+ elseif (strcmpi (str, "maxintervalcount"))
+ maxint = varargin{idx++};
+ elseif (strcmpi (str, "trace"))
+ trace = varargin{idx++};
+ else
+ error ("quadgk: unknown property %s", str);
+ endif
+ else
+ error ("quadgk: expecting property to be a string");
+ endif
+ endwhile
+ if (idx != nargin - 2)
+ error ("quadgk: expecting properties in pairs")
+ endif
+ endif
+ endif
+
+ ## Convert function given as a string to a function handle
+ if (ischar (f))
+ f = @(x) feval (f, x);
+ endif
+
+ ## Use variable subsitution to weaken endpoint singularities and to
+ ## perform integration with endpoints at infinity. No transform for
+ ## contour integrals
+ if (iscomplex (a) || iscomplex (b) || iscomplex(waypoints))
+ ## contour integral, no transform
+ subs = [a; waypoints; b];
+ h = b - a;
+ trans = @(t) t;
+ elseif (isinf (a) && isinf(b))
+ ## Standard Infinite to finite integral transformation.
+ ## \int_{-\infinity_^\infinity f(x) dx = \int_-1^1 f (g(t)) g'(t) dt
+ ## where
+ ## g(t) = t / (1 - t^2)
+ ## g'(t) = 1 / (1 + t^2) ^ 2
+ ## waypoint transform is then
+ ## t = (-1 + sqrt(1 + 4 * g(t) .^ 2)) ./ (2 * g(t))
+ if (!isempty (waypoints))
+ trans = @(x) (-1 + sqrt(1 + 4 * x .^ 2)) ./ (2 * x);
+ subs = [-1; trans(waypoints); 1];
+ else
+ subs = linspace (-1, 1, 11)';
+ endif
+ h = 2;
+ trans = @(t) t ./ (1 - t.^2);
+ f = @(t) f (t ./ (1 + t .^ 2)) ./ ((1 + t .^ 2) .^ 2);
+ elseif (isinf(a))
+ ## Formula defined in Shampine paper as two separate steps. One to
+ ## weaken singularity at finite end, then a second to transform to
+ ## a finite interval. The singularity weakening transform is
+ ## \int_{-\infinity}^b f(x) dx =
+ ## - \int_{-\infinity}^0 f (b - t^2) 2 t dt
+ ## (note minus sign) and the finite interval transform is
+ ## \int_{-\infinity}^0 f(b - t^2) 2 t dt =
+ ## \int_{-1}^0 f (b - g(s) ^ 2 ) 2 g(s) g'(s) ds
+ ## where
+ ## g(s) = s / (1 + s)
+ ## g'(s) = 1 / (1 + s) ^ 2
+ ## waypoint transform is then
+ ## t = sqrt (b - x)
+ ## s = - t / (t + 1)
+ if (!isempty (waypoints))
+ tmp = sqrt (b - waypoints);
+ trans = @(x) - x ./ (x + 1);
+ subs = [0; trans(tmp); 1];
+ else
+ subs = linspace (0, 1, 11)';
+ endif
+ h = 1;
+ trans = @(t) b - (t ./ (1 + t)).^2;
+ f = @(s) - 2 * s .* f (b - (s ./ (1 + s)) .^ 2) ./ ((1 + s) .^ 3);
+ elseif (isinf(b))
+ ## Formula defined in Shampine paper as two separate steps. One to
+ ## weaken singularity at finite end, then a second to transform to
+ ## a finite interval. The singularity weakening transform is
+ ## \int_a^\infinity f(x) dx = \int_0^\infinity f (a + t^2) 2 t dt
+ ## and the finite interval transform is
+ ## \int_0^\infinity f(a + t^2) 2 t dt =
+ ## \int_0^1 f (a + g(s) ^ 2 ) 2 g(s) g'(s) ds
+ ## where
+ ## g(s) = s / (1 - s)
+ ## g'(s) = 1 / (1 - s) ^ 2
+ ## waypoint transform is then
+ ## t = sqrt (x - a)
+ ## s = t / (t + 1)
+ if (!isempty (waypoints))
+ tmp = sqrt (waypoints - a);
+ trans = @(x) x ./ (x + 1);
+ subs = [0; trans(tmp); 1];
+ else
+ subs = linspace (0, 1, 11)';
+ endif
+ h = 1;
+ trans = @(t) a + (t ./ (1 - t)).^2;
+ f = @(s) 2 * s .* f (a + (s ./ (1 - s)) .^ 2) ./ ((1 - s) .^ 3);
+ else
+ ## Davis, Rabinowitz, "Methods of Numerical Integration" p441 2ed.
+ ## Presented in section 5 of the Shampine paper as
+ ## g(t) = ((b - a) / 2) * (t / 2 * (3 - t^2)) + (b + a) / 2
+ ## g'(t) = ((b-a)/4) * (3 - 3t^2);
+ ## waypoint transform can then be found by solving for t with
+ ## Maxima (solve (c + 3*t - 3^3, t);). This gives 3 roots, two of
+ ## which are complex for values between a and b and so can be
+ ## ignored. The third is
+ ## c = (-4*x + 2*(b+a)) / (b-a);
+ ## k = ((sqrt(c^2 - 4) + c)/2)^(1/3);
+ ## t = (sqrt(3)* 1i * (1 - k^2) - (1 + k^2)) / 2 / k;
+ if (! isempty (waypoints))
+ trans = @__quadgk_finite_waypoint__;
+ subs = [-1; trans(waypoints, a, b); 1];
+ else
+ subs = linspace(-1, 1, 11)';
+ endif
+ h = 2;
+ trans = @(t) ((b - a) ./ 4) * t .* (3 - t.^2) + (b + a) ./ 2;
+ f = @(t) f((b - a) ./ 4 .* t .* (3 - t.^2) + (b + a) ./ 2) .* ...
+ 3 .* (b - a) ./ 4 .* (1 - t.^2);
+ endif
+
+ ## Split interval into at least 10 sub-interval with a 15 point
+ ## Guass-Kronrod rule giving a minimum of 150 function evaluations
+ while (length (subs) < 11)
+ subs = [subs' ; subs(1:end-1)' + diff(subs') ./ 2, NaN](:)(1 : end - 1);
+ endwhile
+ subs = [subs(1:end-1), subs(2:end)];
+
+ warn_state = warning ("query", "Octave:divide-by-zero");
+
+ unwind_protect
+ ## Singularity will cause divide by zero warnings
+ warning ("off", "Octave:divide-by-zero")
+
+ ## Initial evaluation of the integrand on the sub-intervals
+ [q_subs, q_errs] = __quadgk_eval__ (f, subs);
+ q0 = sum (q_subs);
+ err0 = sum (q_errs);
+
+ first = true;
+ while (true)
+ ## Check for sub-intervals that are too small. Test must be
+ ## performed in untransformed sub-intervals. What is a good
+ ## value for this test. Shampine suggests 100*eps
+ if (any (diff (trans (subs), [], 2) / (b - a) < 100 * eps))
+ q = q0;
+ err = err0;
+ break;
+ endif
+
+ ## Quit if any evaluations are not finite (Inf or NaN)
+ if (any (! isfinite (q_subs)))
+ warning ("quadgk: non finite integrand encountered");
+ q = q0;
+ err = err0;
+ break;
+ endif
+
+ tol = max (abstol, reltol .* abs (q0));
+
+ ## If the global error estimate is meet exit
+ if (err0 < tol)
+ q = q0;
+ err = err0;
+ break;
+ endif
+
+ ## Accept the sub-intervals that meet the convergence criteria
+ idx = find (abs (q_errs) < tol .* diff (subs, [], 2) ./ h);
+ if (first)
+ q = sum (q_subs (idx));
+ err = sum (q_errs(idx));
+ first = false;
+ else
+ q0 = q + sum (q_subs);
+ err0 = err + sum (q_errs);
+ q += sum (q_subs (idx));
+ err += sum (q_errs(idx));
+ endif
+ subs(idx,:) = [];
+
+ ## If no remaining sub-intervals exit
+ if (rows (subs) == 0)
+ break;
+ endif
+
+ if (trace)
+ disp([fcnt, rows(subs), err, q0]);
+ endif
+
+ ## Split remaining sub-intervals in two
+ mid = (subs(:,2) + subs(:,1)) ./ 2;
+ subs = [subs(:,1), mid; mid, subs(:,2)];
+
+ ## If the maximum sub-interval count is met accept remaining
+ ## sub-interval and exit
+ if (rows (subs) > maxint)
+ warning ("quadgk: maximum interval count (%d) met", maxint);
+ q += sum (q_subs);
+ err += sum (q_errs);
+ break;
+ endif
+
+ ## Evaluation of the integrand on the remaining sub-intervals
+ [q_subs, q_errs] = __quadgk_eval__ (f, subs);
+ endwhile
+
+ if (err > max (abstol, reltol * abs(q)))
+ warning ("quadgk: Error tolerance not met. Estimated error %g", err);
+ endif
+ unwind_protect_cleanup
+ if (strcmp (warn_state.state, "on"))
+ warning ("on", "Octave:divide-by-zero")
+ endif
+ end_unwind_protect
+ endif
+endfunction
+
+function [q, err] = __quadgk_eval__ (f, subs)
+ ## A (15,7) point pair of Guass-Konrod quadrature rules. The abscissa
+ ## and weights are copied directly from dqk15w.f from quadpack
+
+ persistent abscissa = [-0.9914553711208126e+00, -0.9491079123427585e+00, ...
+ -0.8648644233597691e+00, -0.7415311855993944e+00, ...
+ -0.5860872354676911e+00, -0.4058451513773972e+00, ...
+ -0.2077849550078985e+00, 0.0000000000000000e+00, ...
+ 0.2077849550078985e+00, 0.4058451513773972e+00, ...
+ 0.5860872354676911e+00, 0.7415311855993944e+00, ...
+ 0.8648644233597691e+00, 0.9491079123427585e+00, ...
+ 0.9914553711208126e+00];
+
+ persistent weights15 = ...
+ diag ([0.2293532201052922e-01, 0.6309209262997855e-01, ...
+ 0.1047900103222502e+00, 0.1406532597155259e+00, ...
+ 0.1690047266392679e+00, 0.1903505780647854e+00, ...
+ 0.2044329400752989e+00, 0.2094821410847278e+00, ...
+ 0.2044329400752989e+00, 0.1903505780647854e+00, ...
+ 0.1690047266392679e+00, 0.1406532597155259e+00, ...
+ 0.1047900103222502e+00, 0.6309209262997855e-01, ...
+ 0.2293532201052922e-01]);
+
+ persistent weights7 = ...
+ diag ([0.1294849661688697e+00, 0.2797053914892767e+00, ...
+ 0.3818300505051889e+00, 0.4179591836734694e+00, ...
+ 0.3818300505051889e+00, 0.2797053914892767e+00, ...
+ 0.1294849661688697e+00]);
+
+ halfwidth = diff (subs, [], 2) ./ 2;
+ center = sum (subs, 2) ./ 2;;
+ x = bsxfun (@plus, halfwidth * abscissa, center);
+ y = reshape (f (x(:)), size(x));
+
+ ## This is faster than using bsxfun as the * operator can use a
+ ## single BLAS call, rather than rows(sub) calls to the @times
+ ## function.
+ q = sum (y * weights15, 2) .* halfwidth;
+ err = abs (sum (y(:,2:2:end) * weights7, 2) .* halfwidth - q);
+endfunction
+
+function t = __quadgk_finite_waypoint__ (x, a, b)
+ c = (-4 .* x + 2.* (b + a)) ./ (b - a);
+ k = ((sqrt(c .^ 2 - 4) + c) ./ 2) .^ (1/3);
+ t = real ((sqrt(3) .* 1i * (1 - k .^ 2) - (1 + k .^ 2)) ./ 2 ./ k);
+endfunction
+
+%error (quadgk (@sin))
+%error (quadgk (@sin, -pi))
+%error (quadgk (@sin, -pi, pi, 'DummyArg'))
+
+%!assert (quadgk(@sin,-pi,pi), 0, 1e-6)
+%!assert (quadgk(inline('sin'),-pi,pi), 0, 1e-6)
+%!assert (quadgk('sin',-pi,pi), 0, 1e-6)
+%!assert (quadgk(@sin,-pi,pi,'waypoints', 0, 'MaxIntervalCount', 100,
'reltol', 1e-3, 'abstol', 1e-6, 'trace', false), 0, 1e-6)
+%!assert (quadgk(@sin,-pi,pi,1e-6,false), 0, 1e-6)
+
+%!assert (quadgk(@sin,-pi,0), -2, 1e-6)
+%!assert (quadgk(@sin,0,pi), 2, 1e-6)
+%!assert (quadgk(@(x) 1./sqrt(x), 0, 1), 2, 1e-6)
+%!assert (quadgk (@(x) abs (1 - x.^2), 0, 2, 'Waypoints', 1), 2, 1e-6)
+%!assert (quadgk(@(x) 1./(sqrt(x).*(x+1)), 0, Inf), pi, 1e-6)
+%!assert (quadgk (@(z) log (z), 1+1i, 1+1i, 'WayPoints', [1-1i, -1,-1i,
-1+1i]), -pi * 1i, 1e-6)
diff --git a/scripts/general/quadv.m b/scripts/general/quadv.m
new file mode 100644
--- /dev/null
+++ b/scripts/general/quadv.m
@@ -0,0 +1,131 @@
+## Copyright (C) 2008 David Bateman
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} address@hidden =} quadv (@var{f}, @var{a},
@var{b})
+## @deftypefnx {Function File} address@hidden =} quadl (@var{f}, @var{a},
@var{b}, @var{tol})
+## @deftypefnx {Function File} address@hidden =} quadl (@var{f}, @var{a},
@var{b}, @var{tol}, @var{trace})
+## @deftypefnx {Function File} address@hidden =} quadl (@var{f}, @var{a},
@var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{})
+## @deftypefnx {Function File} address@hidden, @var{fcnt}] =} quadl (@dots{})
+##
+## Numerically evaluate integral using adaptive Simpson's rule.
+## @code{quadv (@var{f}, @var{a}, @var{b})} approximates the integral of
+## @address@hidden(@var{x})} to the default absolute tolerance of @code{1e-6}.
+## @var{f} is either a function handle, inline function or string
+## containing the name of the function to evaluate. The function @var{f}
+## must accept a string, and can return a vector representing the
+## approximation to @var{n} different sub-functions.
+##
+## If defined, @var{tol} defines the absolute tolerance to which to
+## which to integrate each sub-interval of @address@hidden(@var{x})}.
+## While if @var{trace} is defined, displays the left end point of the
+## current interval, the interval length, and the partial integral.
+##
+## Additional arguments @var{p1}, etc, are passed directly to @var{f}.
+## To use default values for @var{tol} and @var{trace}, one may pass
+## empty matrices.
+## @seealso{triplequad, dblquad, quad, quadl, quadgk, trapz}
+## @end deftypefn
+
+function [Q, fcnt] = quadv (f, a, b, tol, trace, varargin)
+ if (nargin < 3)
+ print_usage ();
+ endif
+ if (nargin < 4)
+ tol = [];
+ endif
+ if (nargin < 5)
+ trace = [];
+ endif
+ if (isempty (tol))
+ tol = 1e-6;
+ endif
+ if (isempty (trace))
+ trace = 0;
+ endif
+
+ ## Split the interval into 3 abscissa, and apply a 3 point Simpson's rule
+ c = (a + b) / 2;
+ fa = feval (f, a, varargin{:});
+ fc = feval (f, c, varargin{:});
+ fb = feval (f, b, varargin{:});
+ fcnt = 3;
+
+ ## If have edge singularities, move edge point by eps*(b-a) as
+ ## discussed in Shampine paper used to implement quadgk
+ if (isinf (fa))
+ fa = feval (f, a + eps * (b-a), varargin{:});
+ endif
+ if (isinf (fb))
+ fb = feval (f, b - eps * (b-a), varargin{:});
+ endif
+
+ h = (b - a) / 2;
+ Q = (b - a) / 6 * (fa + 4 * fc + fb);
+
+ [Q, fcnt, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, Q, fcnt, abs (b - a),
+ tol, trace, varargin{:});
+
+ if (fcnt > 10000)
+ warning("Maximum iteration count reached");
+ elseif (isnan(Q) || isinf (Q))
+ warning ("Infinite or NaN function evaluations were returned");
+ elseif (hmin < (b - a) * eps)
+ warning ("Minimum step size reached. Possibly singular integral");
+ endif
+endfunction
+
+function [Q, fcnt, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, Q0,
+ fcnt, hmin, tol, trace, varargin)
+ if (fcnt > 10000)
+ Q = Q0;
+ else
+ d = (a + c) / 2;
+ e = (c + b) / 2;
+ fd = feval (f, d, varargin{:});
+ fe = feval (f, e, varargin{:});
+ fcnt += 2;
+ Q1 = (c - a) / 6 * (fa + 4 * fd + fc);
+ Q2 = (b - c) / 6 * (fc + 4 * fe + fb);
+ Q = Q1 + Q2;
+
+ if (abs(a - c) < hmin)
+ hmin = abs (a - c);
+ endif
+
+ if (trace)
+ disp ([fcnt, a, b-a, Q]);
+ endif
+
+ ## Force at least one adpative step
+ if (fcnt == 5 || abs (Q - Q0) > tol)
+ [Q1, fcnt, hmin] = simpsonstp (f, a, c, d, fa, fc, fd, Q1, fcnt, hmin,
+ tol, trace, varargin{:});
+ [Q2, fcnt, hmin] = simpsonstp (f, c, b, e, fc, fb, fe, Q2, fcnt, hmin,
+ tol, trace, varargin{:});
+ Q = Q1 + Q2;
+ endif
+ endif
+endfunction
+
+%!assert (quadv (@sin, 0, 2 * pi), 0, 1e-5)
+%!assert (quadv (@sin, 0, pi), 2, 1e-5)
+
+%% Handles weak singularities at the edge
+%!assert (quadv (@(x) 1 ./ sqrt(x), 0, 1), 2, 1e-5)
+
diff --git a/scripts/general/triplequad.m b/scripts/general/triplequad.m
new file mode 100644
--- /dev/null
+++ b/scripts/general/triplequad.m
@@ -0,0 +1,66 @@
+## Copyright (C) 2008 David Bateman
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} triplequad (@var{f}, @var{xa}, @var{xb},
@var{ya}, @var{yb}, @var{za}, @var{zb}, @var{tol}, @var{quadf}, @dots{})
+## Numerically evaluate a triple integral. The function over with to
+## integrate is defined by @address@hidden, and the interval for the
+## integration is defined by @address@hidden, @var{xb}, @var{ya},
+## @var{yb}, @var{za}, @var{zb}]}. The function @var{f} must accept a
+## vector @var{x} and a scalar @var{y}, and return a vector of the same
+## length as @var{x}.
+##
+## If defined, @var{tol} defines the absolute tolerance to which to
+## which to integrate each sub-integral.
+##
+## Additional arguments, are passed directly to @var{f}. To use the default
+## value for @var{tol} one may pass an empty matrix.
+## @seealso{dblquad, quad, quadv, quadl, quadgk, trapz}
+## @end deftypefn
+
+function Q = triplequad(f, xa, xb, ya, yb, za, zb, tol, quadf, varargin)
+ if (nargin < 7)
+ print_usage ();
+ endif
+ if (nargin < 8 || isempty (tol))
+ tol = 1e-6;
+ endif
+ if (nargin < 9 || isempty (quadf))
+ quadf = @quadgk;
+ endif
+
+ inner = @__triplequad_inner__;
+ if (ischar (f))
+ f = @(x,y,z) feval (f, x, y, z, varargin{:});
+ varargin = {};
+ endif
+
+ Q = dblquad(@(y, z) inner (y, z, f, xa, xb, tol, quadf, varargin{:}),ya, yb,
za, zb, tol);
+endfunction
+
+function Q = __triplequad_inner__ (y, z, f, xa, xb, tol, quadf, varargin)
+ Q = zeros (size(y));
+ for i = 1 : length (y)
+ Q(i) = feval (quadf, @(x) f (x, y(i), z, varargin{:}), xa, xb, tol);
+ endfor
+endfunction
+
+%% These tests are too expensive to run normally. Disable them
+% !#assert (triplequad (@(x,y,z) exp(-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1,
1, [], @quadgk), pi ^ (3/2) * erf(1).^3, 1e-6)
+% !#assert (triplequad (@(x,y,z) exp(-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1,
1, [], @quadl), pi ^ (3/2) * erf(1).^3, 1e-6)
+% !#assert (triplequad (@(x,y,z) exp(-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1,
1, [], @quadv), pi ^ (3/2) * erf(1).^3, 1e-6)
- [Changeset]: Add additional quadrature functions,
David Bateman <=