doubleintegral

Numerical integration of f(x,y) with flexible arguments for giving domain and method.

The intention is to make a more flexible and more general function than the Matlab function dblquad. Especially the domain may be given in a more general way, it may be a circle, a rectangle, or any convex polygon. Also different methods may be used.

Q = DOUBLEINTEGRAL(FUN, DOMAIN, PARAM)
FUN is (normally) a function handle. When method is 'dblquad' the
function Z=FUN(X,Y) should accept a vector X and a scalar Y and return
a vector Z of values as Matlab's 'dblquad'.
For other methods the function just need to handle single element
arguments. The function may accept matrix, and vector, arguments of the
same size and return a matrix also with this size, such that
Z(i) = FUN(X(i),Y(i)). The function will be used this way if parameter
'matrixarg' is also included in the PARAM struct and set to true.
Default is false, and used if 'matrixarg' is not a field in PARAM.
If FUN == 1, then area of domain is returned.

DOMAIN is a struct where the field 'type' gives how the domain is
given. Examples:
domain = struct('type','circle','xc',0,'yc',0,'radius',1);
domain = struct('type','sector','th1',-pi,'th2',pi,'r1',0,'r2',1,'xc',0,'yc',0);
domain = struct('type','box','x1',-1,'x2',1,'y1',-1,'y2',1);
domain = struct('type','area','x1',-1,'x2',1,'y1',@(x) sqrt(x),'y2',@(x) x.^2);
domain = struct('type','polygon','x',[0,2,2,1,0],'y',[0,0,1,2,1]);

PARAM is a struct where the field 'method' always must be included.
It gives the method used, often in both x and y direction.
Examples for Matlab quad, i.e. dblquad, fuction with argument tol,
Gauss Quadrature method and Clenshaw-Curtis method:
param = struct('method','dblquad','tol',1e-6);
param = struct('method','gauss','points',6);
param = struct('method','cc','points',10);
param.matrixarg = true; % use vector-arguments in fun

Q = DOUBLEINTEGRAL(FUN, DOMAIN, PARAM, verbose, makefigure)
verbose = 0 (display only errors), 1 (+ if arguments are overruled), 2
(+ results), 3 (+ some status information), 4 (+ some debugging info)
makefigure = 0 (default) for no figure, 1 for domain, 5 for contour of
function in domain.

Examples:
fun = @(x,y) y*sin(x)+x*cos(y);
domain = struct('type','box','x1',pi,'x2',2*pi,'y1',0,'y2',pi);
param = struct('method','dblquad','tol',1e-6);
Q = doubleintegral(fun, domain, param); % -pi^2

fun = @(x,y) y.*sin(x)+x.*cos(y); % allows vectors in both x and y
domain = struct('type','box','x1',pi,'x2',2*pi,'y1',0,'y2',pi);
param = struct('method','gauss','points',6,'matrixarg',true);
Q = doubleintegral(fun, domain, param, 3, 5); % -pi^2

fun = 1; % Area is one dimensional integral of f(x) = abs(y2(x)-y1(x))
domain = struct('type','area','x1',0,'x2',1,'y1',@(x) x.^2,'y2',@(x) sqrt(x));
param = struct('method','cc','points',12);
Q = doubleintegral(fun, domain, param); % 1/3

fun = @(x,y) sqrt(max(zeros(size(x)), 1-x.^2-y.^2)); % unit sphere
domain = struct('type','circle','xc',0,'yc',0,'radius',1);
param = struct('method','gauss','points',25,'matrixarg',true);
Q = doubleintegral(fun, domain, param); % 2*pi/3
% an alternative is to use dblquad, domain may be a circle or a box.
domain = struct('type','box','x1',-1,'x2',1,'y1',-1,'y2',1);
param = struct('method','dblquad','tol',1e-6);
% or polygon to find just a sector
domain = struct('type','polygon','x',[0,1,1],'y',[0,0,1]);
Q = doubleintegral(fun, domain, param); % pi/12
% or grid of points in sector. Note still fun = @(x,y)
domain = struct('type','sector','th1',0,'th2',pi/4,'r1',0,'r2',1,'xc',0,'yc',0);
param = struct('method','dblquad','tol',1e-6);
Q = doubleintegral(fun, domain, param); % pi/12