function [Q nv] = grnmol( f, V, s ) % % Q = grnmol( f, V ) % computes approximation to the integral of f over an s-simplex % with vertices as columns of V, an n x (n+1) matrix, using % order 1, 2, ..., s (degree 2s+1) Grundmann-Moler rules. % Output Q is a vector approximations of degree 1, 3, ... 2s+1. % Example: % final two results should be 2/(11x10x9x8x7x6) % n = 4; disp( grnmol(@(x)x(1)^2*x(n)^5,eye(n,n+1),4) ) % Reference: % "Invariant Integration Formulas for the N-Simplex by % Combinatorial Methods", A. Grundmann and H. M. Moller, % SIAM J Numer. Anal. 15(1978), pp. 282-290 % [n,np] = size(V); nd = floor(n/2); d = 0; Q = zeros(1,s+1); Qv = Q; Vol = abs(det(V(:,1:n)-V(:,np)*ones(1,n)))/prod([1:n]); while 1, m = n + 2*d + 1; al = ones(n,1); alz = 2*d + 1; Qs = 0; nv = 0; while 1, Qs = Qs + f(V*[alz; al]/m); nv = nv + 1; for j = 1 : n, alz = alz - 2; if alz > 0, al(j) = al(j) + 2; break, end alz = alz + al(j) + 1; al(j) = 1; end, if alz == 2*d+1, break, end end, d = d + 1; Qv(d) = Vol*Qs; Q(d) = 0; p = 2/prod(2*[n+1:m]); for i = 1 : d Q(d) = Q(d) + (m+2-2*i)^(2*d-1)*p*Qv(d+1-i); p = -p*(m+1-i)/i; end, if d > s, break, end end % end grnmol