Mercurial > octave-nkf
changeset 7771:680631e787aa
Add quadv, quadgk, dblquad and triplequad functions
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Sat, 10 May 2008 21:55:40 +0200 |
| parents | c6a1a217ac3c |
| children | ff717f2f9feb |
| files | scripts/ChangeLog scripts/general/Makefile.in scripts/general/dblquad.m scripts/general/quadgk.m scripts/general/quadv.m scripts/general/triplequad.m |
| diffstat | 6 files changed, 707 insertions(+), 4 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog Fri May 09 14:25:46 2008 +0200 +++ b/scripts/ChangeLog Sat May 10 21:55:40 2008 +0200 @@ -1,5 +1,8 @@ 2008-05-12 David Bateman <dbateman@free.fr> + * general/dblquad.m, general/quadgk.m, general/quadv.m, + general/triplequad.m: New functions. + * testfun/assert.m: Allow assert(cond, errmsg, ...) and assert(cond, msgid, errmsg, ...) syntax for compatibility.
--- a/scripts/general/Makefile.in Fri May 09 14:25:46 2008 +0200 +++ b/scripts/general/Makefile.in Sat May 10 21:55:40 2008 +0200 @@ -36,16 +36,16 @@ 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))
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/general/dblquad.m Sat May 10 21:55:40 2008 +0200 @@ -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 @code{@var{f}}, and the interval for the +## integration is defined by @code{[@var{xa}, @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)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/general/quadgk.m Sat May 10 21:55:40 2008 +0200 @@ -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} {[@var{q}, @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 @code{@var{f} (@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 @code{[@var{a}, @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 @code{@var{prop}, @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 @code{@var{x} = 1}. +## +## @item Trace +## If logically true, then @code{quadgk} prints information on the +## convergence of the quadrature at each iteration. +##@end 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)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/general/quadv.m Sat May 10 21:55:40 2008 +0200 @@ -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} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}) +## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}) +## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}) +## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{}) +## @deftypefnx {Function File} {[@var{q}, @var{fcnt}] =} quadl (@dots{}) +## +## Numerically evaluate integral using adaptive Simpson's rule. +## @code{quadv (@var{f}, @var{a}, @var{b})} approximates the integral of +## @code{@var{f}(@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 @code{@var{f}(@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) +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/general/triplequad.m Sat May 10 21:55:40 2008 +0200 @@ -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 @code{@var{f}}, and the interval for the +## integration is defined by @code{[@var{xa}, @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)