anasolv.m

用MATLAB的脚本语言M语言编写的

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;