function Q = gausndv( f, d, n, type, al ) % % Q = gausndv( f, d, n, type, al ) % computes d-dimensional Gauss product rule Q for function f, % using vectorized function calls % Inputs % f integrand function, assuming d-dim column vectors for input % d > 0, number of dimensions % n > 0, number Gauss of points in each dimension % type = 0(default), 1, 2 % Gauss-Legendre, Laguerre, Hermite % al = shape parameter(default = 1) for Laguerre, Hermite % % Example: for integral of exp(-sum(x)^2) over [-1,1]^5, % with 3-point Gauss-Legendre product rule, use % Q = gausndv( @(x)exp(-sum(x).^2), 5, 3 ) % Notes: % a) f needs to allow the input x to be a matrix with d rows, and % f should be vectorized to return a row vector of k values % when x is dxk. % b) if f is defined in an mfile, input name must be quoted 'f'. % Example (same problem as above): % create "gtest.m' mfile (containg the one line) % function f = gtest(x), f = exp(-sum(x).^2); % then use gausndv with the matlab statement % Q = gausndv( 'gtest', 5, 3 ) % or % Q = gausndv(@(x)gtest(x), 5, 3 ) % if nargin < 4, type = 0; end, if nargin < 5, al = 1; end [ x w ] = gausrl( n, type, al ); c = ones(d-1,1); s = zeros(d,1); while 1, s(1) = f([x'; x(c)*ones(1,n)])*w; for i = 1 : d-1, s(i+1) = s(i+1) + w(c(i))*s(i); s(i) = 0; if c(i) == n, c(i) = 1; else, c(i) = c(i) + 1; break, end end, if sum(c) == d-1, break, end end, Q = s(d); % end gaussndv function [ x, w ] = gausrl( n, type, al ), a = zeros(1,n); % % Compute n-point Gauss-rule points x and weights w % Type = 0, Gauss-Legendre, D = [-1,1], w = 1 % Type = 1, Gauss-Laguerre, D = [ 0,inf], w = x^(al-1) exp(-x) % Type = 2, Gauss-Hermite, D = [-inf,inf], w = |x|^(al-1) exp(-x^2) % % Reference: Golub & Welsch, Math Comp 23 (1969), 221-236. % if type == 0, b = 1./sqrt(4-1./[1:n-1].^2); mu = 2; elseif type == 1, if nargin < 3, al = 1; end, mu = gamma(al); a = al + 2*[0:n-1]; b = -sqrt([1:n-1].*([0:n-2]+al)); elseif type == 2, if nargin < 3, al = 1; end b = sqrt([1:n-1]/2); mu = gamma(al/2); if al ~= 1, b(1:2:n-1) = sqrt([0:2:n-2]/2+al/2); end end J = diag(a) + diag(b,1) + diag(b,-1); [ v, d ] = eig(J); x = diag(d); w = mu*v(1,:)'.^2; % % end gausrl %