# dblquad over general domain

### David Woods (view profile)

05 Apr 2006 (Updated )

Extension of dblquad over a rectangle domain to general case.

```function [Q,fcnt] = quadx(funfcn,a,b,tol,trace,varargin)
%   Q = QUAD(FUN,A,B) tries to approximate the integral of function
%   FUN from A to B to within an error of 1.e-6 using recursive
%   adaptive Simpson quadrature.  The function Y = FUN(X) should
%   accept a vector argument X and return a vector result Y, the
%   integrand evaluated at each element of X.
%
%   Q = QUAD(FUN,A,B,TOL) uses an absolute error tolerance of TOL
%   instead of the default, which is 1.e-6.  Larger values of TOL
%   result in fewer function evaluations and faster computation,
%   but less accurate results.  The QUAD function in MATLAB 5.3 used
%   a less reliable algorithm and a default tolerance of 1.e-3.
%
%   [Q,FCNT] = QUAD(...) returns the number of function evaluations.
%
%   QUAD(FUN,A,B,TOL,TRACE) with non-zero TRACE shows the values
%   of [fcnt a b-a Q] during the recursion.
%
%   arguments P1, P2, ... to be passed directly to function FUN,
%   FUN(X,P1,P2,...).  Pass empty matrices for TOL or TRACE to
%   use the default values.
%
%   Use array operators .*, ./ and .^ in the definition of FUN
%   so that it can be evaluated with a vector argument.
%
%   Function QUADL may be more efficient with high accuracies
%   and smooth integrands.
%
%   Example:
%       FUN can be specified as:
%
%       An inline object:
%          F = inline('1./(x.^3-2*x-5)');
%          Q = quad(F,0,2);
%
%       A function handle:
%          Q = quad(@myfun,0,2);
%          where myfun.m is an M-file:
%             function y = myfun(x)
%             y = 1./(x.^3-2*x-5);
%

%   Based on "adaptsim" by Walter Gander.
%   Ref: W. Gander and W. Gautschi, "Adaptive Quadrature Revisited", 1998.
%   http://www.inf.ethz.ch/personal/gander
%   Copyright 1984-2002 The MathWorks, Inc.
%   \$Revision: 5.26 \$  \$Date: 2002/04/08 20:26:50 \$

f = fcnchk(funfcn);
if nargin < 4 | isempty(tol), tol = 1.e-6; end;
if nargin < 5 | isempty(trace), trace = 0; end;

% For the integration over irregular domain,
% the region is likely to be too small, then return 0
% to suppress the warning message. \$Date: 2005/05/07 \$
h=b-a;
if abs(h)<eps
Q=0; fcnt=0;return;
end

% Initialize with three unequal subintervals.
h = 0.13579*h;
x = [a a+h a+2*h (a+b)/2 b-2*h b-h b];
y = feval(f, x, varargin{:});
fcnt = 7;

% Fudge endpoints to avoid infinities.
if ~isfinite(y(1))
y(1) = feval(f,a+eps*(b-a),varargin{:});
fcnt = fcnt+1;
end
if ~isfinite(y(7))
y(7) = feval(f,b-eps*(b-a),varargin{:});
fcnt = fcnt+1;
end

% Call the recursive core integrator.
hmin = eps/1024*abs(b-a);
[Q(1),fcnt,warn(1)] = ...
[Q(2),fcnt,warn(2)] = ...
[Q(3),fcnt,warn(3)] = ...
Q = sum(Q);
warn = max(warn);

switch warn
case 1
'Minimum step size reached; singularity possible.')
case 2
'Maximum function count exceeded; singularity likely.')
case 3
'Infinite or Not-a-Number function value encountered.')
otherwise
% No warning.
end
%fcnt
% ------------------------------------------------------------------------

function [Q,fcnt,warn] = quadstep (f,a,b,fa,fc,fb,tol,trace,fcnt,hmin,varargin)
%QUADSTEP  Recursive core routine for function QUAD.

maxfcnt = 10000;

% Evaluate integrand twice in interior of subinterval [a,b].
h = b - a;
c = (a + b)/2;
if abs(h) < hmin | c == a | c == b
% Minimum step size reached; singularity possible.
Q = h*fc;
warn = 1;
return
end
x = [(a + c)/2, (c + b)/2];
y = feval(f, x, varargin{:});
fcnt = fcnt + 2;
if fcnt > maxfcnt
% Maximum function count exceeded; singularity likely.
Q = h*fc;
warn = 2;
return
end
fd = y(1);
fe = y(2);

% Three point Simpson's rule.
Q1 = (h/6)*(fa + 4*fc + fb);

% Five point double Simpson's rule.
Q2 = (h/12)*(fa + 4*fd + 2*fc + 4*fe + fb);

% One step of Romberg extrapolation.
Q = Q2 + (Q2 - Q1)/15;

if ~isfinite(Q)
% Infinite or Not-a-Number function value encountered.
warn = 3;
return
end
if trace
disp(sprintf('%8.0f %16.10f %18.8e %16.10f',fcnt,a,h,Q))
end

% Check accuracy of integral over this subinterval.
if abs(Q2 - Q) <= tol
warn = 0;
return

% Subdivide into two subintervals.
else