function [hypval,cfs]=hyperlissa(deg,pts,f3); % Stefano De Marchi and Marco Vianello, February 2015 % computes the hyperinterpolation polynomial of a trivariate function % at a rank-1 Chebyshev lattice on a Lissajous curve of the cube % a call has the form [hypval,cfs]=hyperlissa(deg,pts,@f3) where f3.m % is a user defined trivariate Matlab function % NOTE: to be used in connection with the Chebfun package % http://www.chebfun.org % input: % deg: hyperinterpolation degree % pts: 3-column array of evaluation points % output: % hypval: 1-column array of values of the hyperinterpolation polynomial % at the points pts % cfs: hyperinterpolation coefficients n=deg; if mod(n,2)==0 a=3/4*n^2+1/2*n; b=3/4*n^2+n; c=3/4*n^2+3/2*n+1; else a=3/4*n^2+1/4; b=3/4*n^2+3/2*n-1/4; c=3/4*n^2+3/2*n+3/4; end mu=n*c+1; % chebfun along the Lissajous curve g=chebfun(@(t) f3(cos(a*acos(t)),cos(b*acos(t)),cos(c*acos(t))),mu+1); cc=chebpoly(g); % cc=chebcoeffs(g) with Chebfun Version 5 u=pi*cc; u(2:end-1)=u(2:end-1)/2; u=fliplr(u); % indexes i,j,k: graded lexicographical order ind=(0:1:n); [j1,j2,j3]=meshgrid(ind); for s=0:n good=find(j1(:)+j2(:)+j3(:)==s); i(1+s*(s+1)*(s+2)/6:(s+1)*(s+2)*(s+3)/6)=j1(good); j(1+s*(s+1)*(s+2)/6:(s+1)*(s+2)*(s+3)/6)=j2(good); k(1+s*(s+1)*(s+2)/6:(s+1)*(s+2)*(s+3)/6)=j3(good); end; % hyperinterpolation coeffs alpha1=i*a+j*b+k*c+1; alpha2=abs(i*a+j*b-k*c)+1; alpha3=abs(i*a-j*b)+k*c+1; alpha4=abs(abs(i*a-j*b)-k*c)+1; cfs=pi^2/4*sigma(i*a).*sigma(j*b).*sigma(k*c).*(u(alpha1)+u(alpha2)+ ... u(alpha3)+u(alpha4)); % hyperinterpolation values x=pts(:,1);y=pts(:,2);z=pts(:,3); for s=1:length(x) tix=cos(i.*acos(x(s))).*sigma(i); tjy=cos(j.*acos(y(s))).*sigma(j); tkz=cos(k.*acos(z(s))).*sigma(k); val(s)=cfs*(tix.*tjy.*tkz)'; end hypval=val'; end function sigmaval=sigma(m) sigmaval=sqrt(1+sign(m))./sqrt(pi); end