gquad.m

统计分析的软件包

function area = gquad (fun,xlow,xhigh,mparts,bp,wf)
%GQUAD	Legendre Gaussian quadrature.
%	AREA=GQUAD(FUN,XLOW,XHIGH,MPARTS,BP,WF) evaluates the integral
%	of an externally defined function FUN(X) between limits XLOW
%	and XHIGH. The numerical integration is performed using a
%	composite Gauss integration rule.  The whole interval is
%	divided into MPARTS subintervals and the integration over each
%	subinterval is done with an nquad point Gauss formula which
%	involves base points BP and weight factors WF.  The normalized
%	interval of integration for the BP and WF constants is -1 to
%	+1. The algorithm is described by the summation relation
%	x=b                     j=n k=m
%	integral( f(x)*dx ) = d1*sum sum( wf(j)*fun(a1+d*k+d1*bp(j)) )
%	x=a                     j=1 k=1
%	where bp are base points, wf are weight factors, m = mparts,
%	n = length(bp) and d = (b-a)/m, d1 = d/2, a1 = a-d1.  The base
%	points and weight factors must first be generated by a call to
%	grule.  See GRULE.

%	Howard Wilson, U. of Alabama, Spring 1990
%	Revised GKS 5 June 92

bp=bp(:); wf=wf(:);
d = (xhigh - xlow)/mparts;  d2 = d/2;  nquad = length(bp);
x = (d2*bp)*ones(1,mparts) + (d*ones(nquad,1))*(1:mparts);
x = x(:) + (xlow-d2); fv = feval(fun,x); wv = wf*ones(1,mparts);
area=d2*(wv(:)'*fv(:));