derived-bop.zip

function [S,Rts,K]=anasolv(G,cc,dd,x0) G=ss(G); [Ga,xa]=augment(G,cc,dd,x0); vec=ones(1,length(Ga.a)); [Aj,T,K,Rts]=jordan_m(Ga.a,vec); B=inv(T)*xa; F=[]; G0=[]; j=1; for i=1:length(K), if (imag(Rts(j))>1e-5), F1=[B(j),B(j+1); B(j+1),-B(j)]; G1=diag([1,1]); j=j+2; else Xx=B(j:j+K(i)-1)'; F1=Xx; E=1; for k=2:K(i) E(k)=E(k-1)/(k-1); Xx=[Xx(2:K(i)),0]; F1=[F1; Xx]; end G1=diag(E); j=j+K(i); end mm=size(F1); nn=length(F); F=[F,zeros(nn,mm(2));zeros(mm(1),nn),F1]; E=[E,zeros(nn,mm(2));zeros(mm(1),nn),G1]; end S=T*F*G; %jordan_m is a sub function to anasolv() function [AJ,T,JJ,Rt]=jordan_m(A,V0,Eps) if nargin==2, Eps=1e-5; end Rt=eig(A); Rt=sort(Rt); Rt=Rt(length(Rt):-1:1); Nd=length(Rt); M=V0; k=0; JJ=[]; V=[]; for i=2:Nd, M=[M; M(i-1,:)*A]; end while k<Nd, k=k+1; K=0; key=1; Rt0=Rt(k); Rt1=Rt(k); if (abs(imag(Rt0))>Eps), key=2; K=K+1; k=k+2; else while (k<=Nd & abs(Rt0-Rt1)<1e-5), k=k+1; K=K+1; key=0; if k<=Nd, Rt1=Rt(k); end end, end if key~=1, JJ=[JJ K]; k=k-1; end end K=0; for i=1:length(JJ), V0=[]; Rt0=Rt(K+1); j=(1:Nd)'; if (abs(imag(Rt0))<eps),K=K+JJ(i); V0=[V0, real(Rt0).^j]; if JJ(i)>1, V1=V0; for k=2:JJ(i), V2=zeros(Nd,1); for j=k:Nd V2(j)=V1(j-1)*(j-1)/(k-1); end V0=[V0 V2]; V1=V2; end, end else, K=K+2; VV=Rt0.^j; V0=[V0,real(VV),imag(VV)]; end V=[V V0]; end T=inv(M)*V; AJ=inv(T)*A*T; AJ(abs(AJ)<Eps)=0;