scripts/general/dblquad.m

@(
)
@(
)
@(
)
@(
)
@(
)
@(
)
@(
)
author	David Bateman <dbateman@free.fr>
date	Sat, 10 May 2008 21:55:40 +0200
parents
children	cadc73247d65
rev	line source
7771 680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	1 ## Copyright (C) 2008 David Bateman
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	2 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	3 ## This file is part of Octave.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	4 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	5 ## Octave is free software; you can redistribute it and/or modify it
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	6 ## under the terms of the GNU General Public License as published by
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	7 ## the Free Software Foundation; either version 3 of the License, or (at
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	8 ## your option) any later version.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	9 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	10 ## Octave is distributed in the hope that it will be useful, but
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	13 ## General Public License for more details.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	14 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	15 ## You should have received a copy of the GNU General Public License
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	16 ## along with Octave; see the file COPYING. If not, see
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	17 ## <http://www.gnu.org/licenses/>.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	18
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	19 ## -- texinfo --
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	20 ## @deftypefn {Function File} {} dblquad (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{tol}, @var{quadf}, @dots{})
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	21 ## Numerically evaluate a double integral. The function over with to
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	22 ## integrate is defined by @code{@var{f}}, and the interval for the
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	23 ## integration is defined by @code{[@var{xa}, @var{xb}, @var{ya},
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	24 ## @var{yb}]}. The function @var{f} must accept a vector @var{x} and a
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	25 ## scalar @var{y}, and return a vector of the same length as @var{x}.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	26 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	27 ## If defined, @var{tol} defines the absolute tolerance to which to
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	28 ## which to integrate each sub-integral.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	29 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	30 ## Additional arguments, are passed directly to @var{f}. To use the default
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	31 ## value for @var{tol} one may pass an empty matrix.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	32 ## @seealso{triplequad, quad, quadv, quadl, quadgk, trapz}
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	33 ## @end deftypefn
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	34
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	35 function Q = dblquad(f, xa, xb, ya, yb, tol, quadf, varargin)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	36 if (nargin < 5)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	37 print_usage ();
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	38 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	39 if (nargin < 6 \|\| isempty (tol))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	40 tol = 1e-6;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	41 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	42 if (nargin < 7 \|\| isempty (quadf))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	43 quadf = @quadgk;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	44 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	45
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	46 inner = @__dblquad_inner__;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	47 if (ischar (f))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	48 f = @(x,y) feval (f, x, y, varargin{:});
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	49 varargin = {};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	50 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	51
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	52 Q = feval (quadf, @(y) inner (y, f, xa, xb, tol, quadf,
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	53 varargin{:}), ya, yb, tol);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	54 endfunction
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	55
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	56 function Q = __dblquad_inner__ (y, f, xa, xb, tol, quadf, varargin)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	57 Q = zeros (size(y));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	58 for i = 1 : length (y)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	59 Q(i) = feval (quadf, @(x) f(x, y(i), varargin{:}), xa, xb, tol);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	60 endfor
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	61 endfunction
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	62
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	63 %% Nasty integrand to show quadgk off
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	64 %!assert (dblquad (@(x,y) 1 ./ (x+y), 0, 1, 0, 1), 2*log(2), 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	65
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	66 %!assert (dblquad (@(x,y) exp(-x.^2 - y.^2) , -1, 1, -1, 1, [], @quadgk), pi * erf(1).^2, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	67 %!assert (dblquad (@(x,y) exp(-x.^2 - y.^2) , -1, 1, -1, 1, [], @quadl), pi * erf(1).^2, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions David Bateman <dbateman@free.fr> parents: diff changeset	68 %!assert (dblquad (@(x,y) exp(-x.^2 - y.^2) , -1, 1, -1, 1, [], @quadv), pi * erf(1).^2, 1e-6)