Stm.m
上传用户:eighthdate
上传日期:2014-05-24
资源大小:270k
文件大小:1k
源码类别:

其他行业

开发平台:

Matlab

  1. function stm=stm(A);
  2. L = eig(A);
  3. n=length(A);
  4. I=eye(n);
  5. o=ones(n);
  6. H=zeros(n,n);
  7. for i=1:n-1
  8. G(:,i)=L.^i;
  9. end
  10. H=[o(:,1),G];
  11. for i=1:n-1,
  12.   if L(i)==L(i+1),
  13.       for j=1:n,
  14.       H(i+1,j)=(j-1)/L(i+1)*H(i+1,j);
  15.       end
  16.   else
  17.   end
  18. end
  19. K=inv(H);
  20. m=n*n;
  21. for j=1:n
  22. for i=1:n:m
  23. L1(j,i+j-1)=L((i+n-1)/n);
  24. end
  25. end
  26. LI=L1';
  27. FA=zeros(n,m);
  28. for i=1:n:m
  29.    for j=2:n
  30.    FA(:,i:i+n-1)=FA(:,i:i+n-1)+A^(j-1)*K(j,(i+n-1)/n);
  31.    end
  32. F1(:,i:i+n-1)=I*K(1,(i+n-1)/n)+FA(:,i:i+n-1);
  33. end
  34. Head =[
  35. ' The state transition matrix is given by:                              '
  36. '                                                                       '
  37. ' phi(t) = Sum Ci*exp(Li*t) + Sum Dj*t*exp(Lj*t)     i=1,. . .,n-j      '
  38. ' where                                                                 '
  39. ' Li = eigenvalues          & Ci = the corresponding matrix coefficients'
  40. ' Lj = repeated eigenvalues & Dj = the corresponding matrix coefficients'
  41. '                                                                       '];
  42. disp(Head)
  43. %pause(4)
  44. for i=1:n:m
  45. k=(i+n-1)/n;
  46. Li=L(k)
  47.   if k>1
  48.       if L(k)==L(k-1), Di=F1(:,i:i+n-1)
  49.       else  Ci=F1(:,i:i+n-1)
  50.       end
  51.   else  Ci=F1(:,i:i+n-1)
  52.   end
  53. end