nwspgr_Smolyak_smolyak网格_sparsegrid_firstw6q_稀疏网格

%***************************************************************************** %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)<a if (c==d) for i=(c-1):-1:1 c=i; if seq(i)~=0, break, end; end end seq(c) = seq(c)-1; c=c+1; seq(c) = a - sum(seq(1:(c-1))); if (c<d) seq((c+1):d)=zeros(1,d-c); end fs = [fs;seq]; end fs = fs+1; end %************************************************************************************ %SpGrKronProd(): function for generating tensor product quadrature rule %Syntax: % [nodes,weights] = SpGrKronProd(n1d,w1D) %Input: % n1D = cell array of 1D nodes % w1D = cell array of 1D weights %Output: % nodes = matrix of tensor product nodes % weights = vector of tensor product weights %************************************************************************************* function [nodes,weights] = SpGrKronProd(n1D,w1D) nodes = n1D{1} ; weights = w1D{1}; for j=2:length(n1D) newnodes = n1D{j}; nodes = [kron(nodes,ones(size(newnodes,1),1)) kron(ones(size(nodes,1),1),newnodes)]; weights = kron(weights,w1D{j}); end end function [n w] = GQU(l) switch l case 1 n = [5.0000000000000000e-001]; w = [1.0000000000000000e+000]; case 2 n = [7.8867513459481287e-001]; w = [5.0000000000000000e-001]; case 3 n = [5.0000000000000000e-001; 8.8729833462074170e-001]; w = [4.4444444444444570e-001; 2.7777777777777712e-001]; case 4 n = [6.6999052179242813e-001; 9.3056815579702623e-001]; w = [3.2607257743127516e-001; 1.7392742256872484e-001]; case 5 n = [5.0000000000000000e-001; 7.6923465505284150e-001; 9.5308992296933193e-001]; w = [2.8444444444444655e-001; 2.3931433524968501e-001; 1.1846344252809174e-001]; case 6 n = [6.1930959304159849e-001; 8.3060469323313235e-001; 9.6623475710157603e-001]; w = [2.3395696728634746e-001; 1.8038078652407072e-001; 8.5662246189581834e-002]; case 7 n = [5.0000000000000000e-001; 7.0292257568869854e-001; 8.7076559279969723e-001; 9.7455395617137919e-001]; w = [2.0897959183673620e-001; 1.9091502525256090e-001; 1.3985269574463935e-001; 6.4742483084431701e-002]; case 8 n = [5.9171732124782495e-001; 7.6276620495816450e-001; 8.9833323870681348e-001; 9.8014492824876809e-001]; w = [1.8134189168918213