function ps61_6242_2004(n) % function ps61_6242_2004(n) % % solves problem 6.1(b) from 6.242/Fall2004 if nargin<1, n=4; end % default argument tol=1e-13; % tolerance N=10000; % number of samples for H-Infinity norm calculation a=[(1./(1:floor(n/2)))';(2:(ceil(n/2)+1))']; % the vector of 1/k (k=1,...,n) %a=1./(1:n)'; g=(a.^(1/3))./(1+a); % vector of G(1/k) W=1./(repmat(a,1,n)+repmat(a',n,1)); % W [U,D]=eig(W); % W=UDU' where U is orthogonal, D diagonal d=diag(D); d=max(d,tol); % regularization ug=U'*g; % solving for c c=U*(ug./d); H2Ghat=sum((abs(ug).^2)./d); % square of H2 norm of Ghat %c=W\g; H2Ghat=g'*c; w=linspace(0,10,N)'; % column vector of frequency samples Ghat=(1./(repmat(a',N,1)+repmat(j*w,1,n)))*c; % samples of Ghat G=((j*w).^(1/3))./(1+j*w); % samples of G disp(['numeric H2 error: ',num2str((10/N)*sum(abs(G-Ghat).^2))]) disp(['H2 error: ',num2str(1-H2Ghat)]) disp(['H-Infinity error: ',num2str(max(abs(G-Ghat)))]) close(gcf) subplot(2,1,1); plot(w,real(G),w,real(Ghat)); grid subplot(2,1,2); plot(w,imag(G),w,imag(Ghat)); grid