function pvldemo(action) %PVLDEMO Padé via Lanczos (PVL) Demostration % Zhaojun Bai % Professor of Computer Science and Mathematics % 3005 Kemper Hall % Department of Computer Science % University of California % One Shields Avenue % Davis, CA 95616 % % Phone: (530) 752-4874 % Fax: (530) 752-4767 % email: bai@cs.ucdavis.edu % % $Revision: 1.10 $ $Date: 2004/03/19 18:30:30 $ if nargin<1, action='initialize'; end; % Problem Dependent Parameters bf_awebook = 'AWEBOOK is taken from ..., ... is shown in the above axes, ... and ... '; ep_awebook = 'Expansion Point: pi*1e9'; or_awebook = ' Order 4| Order 5| Order 6'; mat_awebook = str2mat('4','5','6'); bf_brake834 = 'BRAKE384 is taken from ..., ... is shown in the above axes, ... and ... '; ep_brake834 = 'Expansion Point: 0'; or_brake834 = ' Order 39| Order 40| Order 41'; mat_brake834 = str2mat('39','40','41'); bf_peec = 'PEEC is taken from ..., ... is shown in the above axes, ... and ... '; ep_peec = 'Expansion Point: 2*pi*1e9'; or_peec = ' Order 50| Order 55| Order 60| Order 65 | Order 70'; mat_peec = str2mat('50','55','60','65','70'); if strcmp(action,'initialize'), oldFigNumber=watchon; figNumber=figure( ... 'Visible','off', ... 'NumberTitle','off', ... 'MenuBar', 'none', ... 'Name','Padé via Lanczos (PVL) Demo'); pos=get(figNumber,'Position'); width=pos(3)-pos(1); pos(1)=pos(1)-width*0.5; pos(3)=pos(3)+width*0.5; set(figNumber,'Position',pos); %=================================== % Set up the axes left=0.09; bottom=0.10; width=0.60; height=0.38; axRel=axes( ... 'Units','normalized', ... 'Position',[left bottom width height]); axMag=axes( ... 'Units','normalized', ... 'Position',[left 2*bottom+height width height]); %==================================== % Information for all buttons labelColor=[0.8 0.8 0.8]; top=0.95; bottom=0.05; left=0.75; yInitLabelPos=0.90; left=0.75; labelWid=0.20; labelHt=0.05; btnWid=0.20; btnHt=0.05; % Spacing between the label and the button for the same command btnOffset=0.003; % Spacing between the button and the next command's label spacing=0.05; %==================================== % The CONSOLE frame frmBorder=0.02; yPos=0.05-frmBorder; frmPos=[left-frmBorder yPos btnWid+2*frmBorder 0.9+2*frmBorder]; uicontrol( ... 'Style','frame', ... 'Units','normalized', ... 'Position',frmPos, ... 'BackgroundColor',[0.50 0.50 0.50]); %==================================== % The EXAMPLE popup button btnNumber=1; yLabelPos=top-(btnNumber-1)*(btnHt+labelHt+spacing); labelStr='Example'; labelList=' awebook| brake834| peec'; cmdList=str2mat( ... 'awebook','brake834','peec'); callbackStr='pvldemo(''reset'')'; % Generic label information labelPos=[left yLabelPos-labelHt labelWid labelHt]; uicontrol( ... 'Style','text', ... 'Units','normalized', ... 'Position',labelPos, ... 'BackgroundColor',labelColor, ... 'HorizontalAlignment','left', ... 'String',labelStr); % Generic popup button information btnPos=[left yLabelPos-labelHt-btnHt-btnOffset btnWid btnHt]; hndlEx=uicontrol( ... 'Style','popup', ... 'Units','normalized', ... 'Position',btnPos, ... 'String',labelList, ... 'Callback',callbackStr, ... 'UserData',cmdList); %==================================== % The ORDER popup button btnNumber=2; yLabelPos=top-(btnNumber-1)*(btnHt+labelHt+spacing); labelStr='Padé approximation'; labelList=or_awebook; cmdList=mat_awebook; callbackStr='pvldemo(''eval'')'; % Generic label information labelPos=[left yLabelPos-labelHt labelWid labelHt]; uicontrol( ... 'Style','text', ... 'Units','normalized', ... 'Position',labelPos, ... 'BackgroundColor',labelColor, ... 'HorizontalAlignment','left', ... 'String',labelStr); % Generic popup button information btnPos=[left yLabelPos-labelHt-btnHt-btnOffset btnWid btnHt]; hndlOr=uicontrol( ... 'Style','popup', ... 'Units','normalized', ... 'Position',btnPos, ... 'String',labelList, ... 'Callback',callbackStr, ... 'UserData',cmdList); %==================================== % The EXPANSION POINT text field btnNumber=3; yLabelPos=top-(btnNumber-1)*(btnHt+labelHt+spacing); labelStr=ep_awebook; % Generic label information labelPos=[left yLabelPos-labelHt labelWid labelHt]; hndlEp=uicontrol( ... 'Style','text', ... 'Units','normalized', ... 'Position',labelPos, ... 'BackgroundColor',labelColor, ... 'HorizontalAlignment','left', ... 'String',labelStr); %==================================== % The static text field btnNumber=4; yLabelPos=top-(btnNumber-1)*(btnHt+labelHt+spacing); labelColor=[0.94 0.94 0.94]; labelStr=bf_awebook; % Generic label information labelPos=[left yLabelPos-labelHt-2*btnHt labelWid labelHt+3*btnHt]; hndlBf=uicontrol( ... 'Style','text', ... 'Units','normalized', ... 'Position',labelPos, ... 'BackgroundColor',labelColor, ... 'HorizontalAlignment','left', ... 'String',labelStr); %==================================== % The info button uicontrol( ... 'Style','pushbutton', ... 'Units','normalized', ... 'Position',[left bottom+2*btnHt+spacing btnWid 2*btnHt], ... 'String','Info', ... 'Callback','pvldemo(''info'')'); % The close button. uicontrol( ... 'Style','pushbutton', ... 'Units','normalized', ... 'Position',[left bottom btnWid 2*btnHt], ... 'String','Close', ... 'Callback','close(gcf)'); % Uncover the figure hndlList=[axRel axMag hndlEx hndlOr hndlEp hndlBf]; watchoff(oldFigNumber); set(figNumber,'Visible','on', ... 'UserData',hndlList); pvldemo('eval'); elseif strcmp(action,'info'), helpwin(mfilename) elseif strcmp(action,'reset'), hndlList=get(gcf,'UserData'); hndlEx = hndlList(3); hndlOr = hndlList(4); hndlEp = hndlList(5); hndlBf = hndlList(6); StrList = get(hndlEx, 'UserData'); StrVal = get(hndlEx, 'Value'); Str = deblank(StrList(StrVal,:)); if strcmp(Str, 'awebook'), set(hndlOr, 'String', or_awebook, 'UserData', mat_awebook); set(hndlEp, 'String', ep_awebook); set(hndlBf, 'String', bf_awebook); elseif strcmp(Str, 'brake834'), set(hndlOr, 'String', or_brake834, 'UserData', mat_brake834); set(hndlEp, 'String', ep_brake834); set(hndlBf, 'String', bf_brake834); elseif strcmp(Str, 'peec'), set(hndlOr, 'String', or_peec, 'UserData', mat_peec); set(hndlEp, 'String', ep_peec); set(hndlBf, 'String', bf_peec); end pvldemo('eval'); elseif strcmp(action,'eval'), % Assemble and execute the completed command hndlList=get(gcf,'UserData'); % order of reduced model hndlOr = hndlList(4); StrList = get(hndlOr, 'UserData'); StrVal = get(hndlOr, 'Value'); Str = deblank(StrList(StrVal,:)); n = str2num(Str); hndlEx = hndlList(3); StrList = get(hndlEx, 'UserData'); StrVal = get(hndlEx, 'Value'); Str = deblank(StrList(StrVal,:)); if strcmp(Str, 'awebook'), %==================================== % AWEBOOK % read in data awebookeg; N = size(C,1); % expansion point s0 = 1*pi*1e+9; % Reduced-order modeling [Cn,Gn,bn,ln] = romlin(N,C,G,b,l,n,s0); % Bode plot npts = 1001; freqlow = 10^0; freqhig = 2*10^(10); linst = (freqhig-freqlow)/(npts-1); freq = freqlow:linst:freqhig; s = 2*pi*sqrt(-1)*freq; freqpts = max(size(s)); for k = 1:freqpts, Hn(k) = ln'*((Gn + s(k)*Cn)\bn); end % Transfer function - exact for k = 1:freqpts, HN(k) = l'*((G+s(k)*C)\b); end % Transfer function - approximation for k = 1:freqpts, Hn(k) = ln'*((Gn + s(k)*Cn)\bn); end axMag = hndlList(2); axes(axMag); plot(abs(freq),20*log10(abs(HN)),'r-', abs(freq),20*log10(abs(Hn)),'b--'); xlabel('Frequency (Hz)', 'fontsize', 12); ylabel('Gain (dB)', 'fontsize', 12); legend('Exact', 'PVL'); %set(gca, 'fontsize', 12, 'linewidth', 1.5); axRel = hndlList(1); axes(axRel); plot(abs(freq),angle(HN)*180/pi,'r-', abs(freq),angle(Hn)*180/pi,'b--'); xlabel('Frequency (Hz)', 'fontsize', 12); ylabel('phase(degree)', 'fontsize', 12); legend('Exact', 'PVL'); %set(gca, 'fontsize', 12, 'linewidth', 1.5); elseif strcmp(Str, 'brake834'), % load matrix data, and exact values of transfer function load brake834M; load brake834K; C = M; G = K; N = size(C,1); l = zeros(N,1); l(1)=1; b = l; % expansion point s0 = 0; % Reduced-order modeling [Cn,Gn,bn,ln] = romlin(N,C,G,b,l,n,s0); % Bode plot npts = 501; freqlow = 0; freqhig = 10000; freq = linspace(freqlow,freqhig,npts); s = 2*pi*sqrt(-1)*freq; freqpts = max(size(s)); % evaluate the transfer function of the reduced-order model for k = 1:freqpts, Hn(k) = ln'*((Gn+s(k)^2*Cn)\bn); end % ``exact'' transfer function %for k = 1:freqpts, % HN(k) = l'*((G+s(k)^2*C)\b); %end % precalculated load brake834freqexact; axMag = hndlList(2); axes(axMag); semilogy(abs(freq),abs(HN),'r-', abs(freq),abs(Hn),'b--'); xlabel('Frequency (Hz)', 'fontsize', 12); ylabel('Magnitude', 'fontsize', 12); legend('Exact', 'PVL'); %set(gca, 'fontsize', 12, 'linewidth', 1.5); axRel = hndlList(1); axes(axRel); semilogy(abs(freq),abs(HN-Hn)./abs(HN),'b-'); xlabel('Frequency (Hz)', 'fontsize', 12); ylabel('Relative error', 'fontsize', 12); %set(gca, 'fontsize', 12, 'linewidth', 1.5); elseif strcmp(Str, 'peec'), % load data load ex306v4; % matrix data, and exact values of transfer function C = CM; G = GM; l = LM(:,1); b = LM(:,2); N=size(C,1); % expansion point s0 = 2*pi*1e+9; % reduced-order modeling [Ck,Gk,bk,lk]=romlin(N,C,G,b,l,n,s0); % Bode plot npts = 1001; freqlow = 10^0; freqhig = 0.5*10^(10); % linear space linst = (freqhig-freqlow)/(npts-1); freq = freqlow:linst:freqhig; s = 2*pi*sqrt(-1)*freq; freqpts = max(size(s)); % ``exact'' transfer function % tic; % for k = 1:freqpts, % HN(k) = l'*((G+s(k)*C)\b); % end % Fullmodel_cputime=toc % precalculated load ex306v4HN; % transfer function of reduced-order model tic; for k = 1:freqpts, Hn(k) = lk'*((Gk+s(k)*Ck)\bk); end Reducedordermodel_cputime = toc; % This is not exactly bode polt, but commonly as presented in literature axMag = hndlList(2); axes(axMag); plot(abs(freq),abs(HN),'r-', abs(freq),abs(Hn),'b--'); xlabel('Frequency (Hz)', 'fontsize', 12); ylabel('Voltage (Volts)', 'fontsize', 12); legend('Exact', 'PVL'); %set(gca, 'fontsize', 12, 'linewidth', 1.5); axRel = hndlList(1); axes(axRel); semilogy(abs(freq),abs(HN-Hn)./abs(HN),'b-'); xlabel('Frequency (Hz)', 'fontsize', 12); ylabel('Relative error', 'fontsize', 12); %set(gca, 'fontsize', 12, 'linewidth', 1.5); end end; % if strcmp(action, ...