%***************************************************************************** %nwSpGr: Nodes and weights for numerical integration on sparse grids (Smolyak) %(c) 2007 Florian Heiss, Viktor Winschel %Nov 11, 2007 %***************************************************************************** %************************************************************************************ %nwspgr(): main function for generating nodes & weights for sparse grids intergration %Syntax: % [n w] = nwSpGr(type, dim, k, sym) %Input: % type = String for type of 1D integration rule: % "KPU": Nested rule for unweighted integral over [0,1] % "KPN": Nested rule for integral with Gaussian weight % "GQU": Gaussian quadrature for unweighted integral over [0,1] (Gauss-Legendre) % "GQN": Gaussian quadrature for integral with Gaussian weight (Gauss-Hermite) % func: any function name. Function must accept level l and % return nodes n and weights w for univariate quadrature rule with % polynomial exactness 2l-1 as [n w] = feval(func,level) % Example: If you have installed the CompEcon Toolbox of Paul Fackler and Mario Miranda % (http://www4.ncsu.edu/~pfackler/compecon/toolbox.html) % and you want to integrate with a lognormal weight function, you can use % nwSpGr('qnwlogn', dim, k) % dim : dimension of the integration problem % k : Accuracy level. The rule will be exact for polynomial up to total order 2k-1 % sym : (optional) only used for own 1D quadrature rule (type not "KPU",...). If % sym is supplied and not=0, the code will run faster but will % produce incorrect results if 1D quadrature rule is asymmetric %Output: % n = matrix of nodes with dim columns % w = row vector of corresponding weights %************************************************************************************* function [nodes, weights] = nwspgr(type, dim, k, sym) % interpret inputs if (nargin < 3); error('not enough input arguments') end builtinfct = (strcmp(type,'GQU') || strcmp(type,'GQN') || strcmp(type,'KPU') || strcmp(type,'KPN')); if (nargin == 3) if builtinfct; sym = true; else sym = false; end else sym = logical(sym); end % get 1D nodes & weights try n1D = cell(k,1); w1D = cell(k,1); R1D = zeros(1,k); for level=1:k [n w] = feval(type,level); % if user supplied symmetric 1D rule: keep only positive orthant if ~builtinfct && sym numnew = rows(n); [n sortvec] = sortrows(n); w = w(sortvec); n = n((floor(numnew/2)+1):numnew,:); w = w((floor(numnew/2)+1):numnew,:); end R1D(level) = length(w); n1D{level} = n; w1D{level} = w; end catch error('Error evaluating the 1D rule'); end % initialization minq = max(0,k-dim); maxq = k-1; nodes = []; weights = []; % outer loop over q for q = minq:maxq r = length(weights); bq = (-1)^(maxq-q) * nchoosek(dim-1,dim+q-k); % matrix of all rowvectors in N^D_{q} is = SpGrGetSeq(dim,dim+q); % preallocate new rows for nodes & weights Rq = prod(R1D(is),2); sRq = sum(Rq); nodes = [nodes ; zeros(sRq,dim)]; weights = [weights ; zeros(sRq,1)]; % inner loop collecting product rules for j=1:size(is,1) midx = is(j,:); [newn,neww] = SpGrKronProd(n1D(midx),w1D(midx)); nodes((r+1):(r+Rq(j)),:) = newn; weights((r+1):(r+Rq(j))) = bq .* neww; r = r + Rq(j); end % collect equal nodes: first sort [nodes sortvec] = sortrows(nodes); weights = weights(sortvec); keep = 1; lastkeep = 1; % then make a list of rows to keep and sum weights of equal nodes for j=2:size(nodes,1) if nodes(j,:)==nodes(j-1,:) weights(lastkeep)=weights(lastkeep)+weights(j); else lastkeep = j; keep = [keep ; j ]; end end % drop equal rows nodes = nodes(keep,:); weights = weights(keep); end % If symmetric rule: so far, it's for the positive orthant. Copy to other % orthants! if sym nr = length(weights); m = n1D{1}; for j=1:dim keep = zeros(nr,1); numnew = 0; for r=1:nr if nodes(r,j) ~= m numnew=numnew+1; keep(numnew)=r; end end if numnew>0 nodes = [nodes ; nodes(keep(1:numnew),:)]; nodes(nr+1:nr+numnew,j) = 2*m - nodes(nr+1:nr+numnew,j); weights = [weights ; weights(keep(1:numnew))]; nr=nr+numnew; end end % 3. final sorting [nodes sortvec] = sortrows(nodes); weights = weights(sortvec); end % normalize weights to account for rounding errors weights = weights / sum(weights); end %************************************************************************************** %SpGrGetSeq(): function for generating matrix of all rowvectors in N^D_{norm} %Syntax: % out = nwSpGr(d,norm) %Input: % d = dimension, will be #columns in output % norm = row sum of elements each of the rows has to have %Output: % out = matrix with d columns. Each row represents one vector with all elements >=1 % and the sum of elements == norm %************************************************************************************** function fs = SpGrGetSeq(d, norm) seq = zeros(1,d); a=norm-d; seq(1)=a; fs = seq; c=1; while seq(d)