multiquad.m

工具箱 (数值积分工具箱) matlab 中使用

function [Q,E] = multiquad(funfcn,s,order)
%MULTIQUAD   Multiple integration using a lattice method.
%   Q = MULTIQUAD(F,S) approximates the S-dimensional integral of a  
%   function F(X) over the unit cube, with 3 <= S <= 8.
%   F(X) returns a column vector of values for a matrix-valued input
%   (where each row of the input matrix is a point in S-dimensional space).
%
%   [Q,E] = MULTIQUAD(F,S) returns an estimate of the error.
%
%   [Q,E] = MULTIQUAD(F,S,ORDER)  uses different numbers of function evaluations:
%     ORDER=1 uses approximately   5000 evaluations (this is the default),
%     ORDER=2 uses approximately  20000 evaluations,
%     ORDER=3 uses approximately  80000 evaluations, and
%     ORDER=4 uses approximately 320000 evaluations.
%
% For example,
%  >> f=inline('prod(1+2*pi^2*(x.^2-x+1/6),2)','x');
%  >> [Q,E]=multiquad(f,6,2)
% returns Q = 0.9513 and E = 0.1063. The exact value is Q=1.

% Author: Robert Piche, Tampere University of Technology, 25.4.2003
% Ref: Stephen Joe & Ian H. Sloan, ACM Trans. on Math. Soft. 
%    vol 19 No 4 1993 p. 523-545 and vol 20 No 2 1994 p. 245.

if nargin < 3, order=1; end
if ~ismember(order,1:4), error('invalid ORDER'); end
if (s<3|s>8), error('invalid number of dimensions S'); end

% method constants
n=2;
zList1={NaN,NaN,[441,1,115],[51,1,97,252],[1,116,61,110,11],...
   [41,10,15,33,22,1],[1,25,5,13,2,24,38],[16,14,7,4,8,13,2,1]};
zList2={NaN,NaN,[1,1868,1025],[932,732,1,569],[216,180,1,125,150],...
        [223,290,1,275,248,192],[106,109,57,90,93,76,1],...
		[1,17,27,18,12,8,65,58]};
zList3={NaN,NaN,[2325,1,1845],[1889,1,792,191],[1514,1951,792,151,1],...
       [1,175,649,440,1165,288],[542,289,161,1,71,358,602],...
	   [85,83,205,70,1,265,176,113]};
zList4={NaN,NaN,[28219,5361,1],[15513,1,1158,7318],[4443,9179,1,4331,6445],...
       [4718,4538,1,1229,162,3981],[728,2289,1,2097,1897,1851,914],...
	   [294,903,30,606,806,255,155,1]};
zList={zList1,zList2,zList3,zList4};

mList=[NaN,NaN,619,313,157,79,41,19
       NaN,NaN,2503,1249,619,313,157,79
       NaN,NaN,10007,5003,2503,1249,619,313
       NaN,NaN,39989,19997,10007,5003,2503,1249];

% initialize some variables
z=zList{order}{s};
m=mList(order,s);
k=zeros(1,s);
Q=zeros(1,s+1);
Qu=zeros(1,s);
k(1)=-1;
l=1;
omega=0;

% main body of algorithm
while l<=s
	if k(l)<n-1
		k(l)=k(l)+1;
		l=1;
		val=0;
		j=(0:m-1)';
		x=mod(j*z/m+repmat(k,m,1)/n,1);
		% apply a periodizing transformation to f
		f=prod(6*x.*(1-x),2).*feval(funfcn,(3-2*x).*(x.^2)); 
		val=sum(f); 
		Q(omega+1:s+1)=Q(omega+1:s+1)+val;
		ii=find(k==0);
		Qu(ii)=Qu(ii)+val;
		if omega==0, omega=1; end
	else
		k(l)=0;
		l=l+1;
		if omega<l, omega=l; end
	end
end
Q=(Q./(n.^(0:s)))/m;
Q=Q(s+1);
Qu=Qu/(n^(s-1))/m;
E=norm(Q-Qu)/sqrt(s);