function [w x]=squaresparse1dtomd(d,r,rule) ## Converts a single variable rule to a multidimensional rule on a square, cube, etc ## [w x]=squaresparse1dtomd(d,r,rule) ## ## Inputs ## d: number of variables for the final rule ## r: the amount of refinement ## rule: the class of rules to use options are: 'Clenshaw Curtis' or 'Gauss Legendre' ## ## Outputs ## w: The weights of the sparse grid multidimensional rule ## x: The nodes of the sparse grid multidimensional rule ##r=3; ##d=3; ##rule='gauss legendre';%'clenshaw curtis';%'clenshaw curtis';% rule=lower(rule); wx1d=struct('nodes',[],'weights',[]); for i=1:r+1 switch rule case 'clenshaw curtis' if i==1 [weights,nodes]=clenshawcurtiswx(i); else [weights,nodes]=clenshawcurtiswx(2^(i-1)+1); end case 'gauss legendre' [weights,nodes]=gausslegendrewx(i); endswitch wx1d(i).nodes=nodes; wx1d(i).weights=weights; end w=[]; x=[]; for l = max(d,r+1):d+r if l == d x=repmat(wx1d(1).nodes,1,d); w=prod(repmat(wx1d(1).weights,1,d),2); w=(-1)^(d+r-l)*nchoosek(d-1,d+r-l)*w; else alpha=ones(1,d); alpha(1)=l-d+1; count=nchoosek(l-1,d-1); for i=1:count tmpnodes=wx1d(alpha(1)).nodes; tmpweights=wx1d(alpha(1)).weights; for j=2:d tmpnodes2=wx1d(alpha(j)).nodes; tmpweights2=wx1d(alpha(j)).weights; tmpnodes=[kron(ones(size(tmpnodes2,1),1),tmpnodes) kron(tmpnodes2,ones(size(tmpnodes,1),1))]; tmpweights=kron(tmpweights2,tmpweights); end tmpweights=(-1)^(d+r-l)*nchoosek(d-1,d+r-l)*tmpweights; x=[x;tmpnodes]; w=[w;tmpweights]; if i=3 scatter3(x(:,1),x(:,2),x(:,3),[],w); axis equal; elseif d==2 scatter(x(:,1),x(:,2),[],w); axis equal; elseif d==1 plot(x,w,'o') axis equal; end end function res=nextrestrictedcomposition(c); ## Returns the next restricted composition in lexicographical order. i=1; while c(i)==1 i=i+1; end t=c(i); c(i)=1; c(1)=t-1; c(i+1)=c(i+1)+1; res=c; end function [w,x]=uniquewx(w,x,tol) ## Removes duplicated x and sums corisponding w. points=rows(w); d=columns(x); Iall=zeros(points,1); I=zeros(0,1); J=zeros(0,1); i=1; j=1; ii=1; ww=zeros(0,1); for i=1:points if any(Iall==i) continue else eq=sumsq(x-x(i,:),2).^.5