function [pts,w,momerr] = comprexcub(deg,X,omega,pos) % compression of bivariate cubature formulas by Tchakaloff points % or approximate Fekete points % useful, for example, in node reduction of algebraic cubature formulas % see the web page: http://www.math.unipd.it/~marcov/Tchakaloff.html % by Federico Piazzon, Alvise Sommariva and Marco Vianello % University of Padova, May 2016 % INPUT: % deg: polynomial exactness degree % X: 2-column array of point coordinates % omega: 1-column array of weights % pos: NNLS for pos=1, QR with column pivoting for pos=0 % OUTPUT: % pts: 2-column array of extracted points % w: 1-column array of corresponding weights (positive for pos=1) % momerr: moment reconstruction error % FUNCTION BODY % total-degree Chebyshev-Vandermonde matrix at X rect=[min(X(:,1)) max(X(:,1)) min(X(:,2)) max(X(:,2))]; V=chebvand(deg,X,rect); % polynomial basis orthogonalization [Q,R]=qr(V,0); % tiny imaginary parts could appear Q=real(Q); % possible re-orthogonalization % [Q,R]=qr(Q,0); % moments of the orthogonal basis orthmom=Q'*omega; % weigths computation if pos == 1 % Tchakaloff points (positive weights) weights=lsqnonneg(Q',orthmom); else % approximate Fekete points (possible negative weights) weights=Q'\orthmom; end; % indexes of nonvanishing weights and compression ind=find(abs(weights)>0); pts=X(ind,:); w=weights(ind); % moment reconstruction error % bivariate Chebyshev basis % mom=V'*omega; % momerr=norm(V(ind,:)'*w-mom); % discrete OP basis momerr=norm(Q(ind,:)'*w-orthmom); end function V = chebvand(deg,X,rect); % INPUT: % deg = polynomial degree % X = 2-column array of point coordinates % rect = 4-component vector such that the rectangle % [rect(1),rect(2)] x [rect(3),rect(4)] contains X % OUTPUT: % V = Chebyshev-Vandermonde matrix at X, graded lexic. order % FUNCTION BODY % rectangle containing the mesh if isempty(rect) rect=[min(X(:,1)) max(X(:,1)) min(X(:,2)) max(X(:,2))]; end; % couples with length less or equal to deg % graded lexicographical order j=(0:1:deg); [j1,j2]=meshgrid(j); dim=(deg+1)*(deg+2)/2; couples=zeros(dim,2); for s=0:deg good=find(j1(:)+j2(:)==s); couples(1+s*(s+1)/2:(s+1)*(s+2)/2,:)=[j1(good) j2(good)]; end; % mapping the mesh in the square [-1,1]^2 a=rect(1);b=rect(2);c=rect(3);d=rect(4); map=[(2*X(:,1)-b-a)/(b-a) (2*X(:,2)-d-c)/(d-c)]; % Chebyshev-Vandermonde matrix on the mesh T1=chebpolys(deg,map(:,1)); T2=chebpolys(deg,map(:,2)); V=T1(:,couples(:,1)+1).*T2(:,couples(:,2)+1); end function T=chebpolys(deg,x) % computes the Chebyshev-Vandermonde matrix on the real line by recurrence % INPUT: % deg = maximum polynomial degree % x = 1-column array of abscissas % OUTPUT % T: Chebyshev-Vandermonde matrix at x, T(i,j+1)=T_j(x_i), j=0,...,deg T=zeros(length(x),deg+1); t0=ones(length(x),1); T(:,1)=t0; t1=x; T(:,2)=t1; for j=2:deg t2=2*x.*t1-t0; T(:,j+1)=t2; t0=t1; t1=t2; end; end