function ps7_3_2(n,k) %Attempt to find a reduced model with stable part of %order k for a Riemann sum approximation of order n. C=ones(1,100); %Create a state-space description of the Ga(s) B=[]; a=[]; for m=0:99 B=[B; 1/(200+m)]; a=[a -1*(1+m/100)]; end A=diag(a); D=0; Ga=pck(A,B,C,D); [sysout,sig]=sysbal(Ga,1e-20); %Balance the system. The 1e-20 can be set to %any appropriate tolerance param that gives %enough poles for hankmr to work with! if (k>length(sig)) k=length(sig)-1; %can't choose k greater than what sysbal will end %give us! sh=hankmr(sysout,sig,k,'d'); %model reduce! [Ar,Br,Cr,Dr]=unpck(sh); Ghat=tf(ss(Ar,Br,Cr,Dr)); w=logspace(-4,4,200); %numerically compute the infinity norm of %G-Ghat. Ghattemp=freqresp(Ghat,w) ; %Frequency response of Ghat. for m=1:200 Ghatfreq(m)=Ghattemp(:,:,m); %dealing with strange MATLAB notation... end o=ones(1,200); %Freq response of G minus unstable term. Gfreq=o./(j*w-o).*log(2*o.*(j*w+2*o)./(3*o)./(j*w+o)); Gdif=Gfreq+Ghatfreq; %the plus is correct---remember the minus %sign from the partial fraction expansion! semilogx(w,abs(Gdif)); %frequency response of difference grid title('|G(jw)-Ghat(jw)|') %plot the frequency response to check the xlabel('w') %infinity norm of the difference. s=tf('s') Ghattf=log(1.5)/(s-1)-Ghat %the final transfer function