The ABCs of RBCs: matlab程序问题，达人请指教

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

第六章136页，脚注，下载一个一般化的Schur decomposition 的程序。下载下来后在 7.6版的matlab里面运行出错，说找不到程序。路径没问题，估计是版本问题，anderson这个程序很老了，下载下来的是一个m文件和dll文件。http://www.math.niu.edu/~anderson/。 不知道哪位高人用过这个程序，指点一下!


此外我试图自己用 matlab内置的 ordqz来分解，但根本读不懂使用说明，郁闷！

具体问题很简单，就是将 A，B两个矩阵 进行QZ分解，但要求特征值是按照其绝对值的升序排列。qz直接分解的结果不符合要求。但ordqz不会用。

B=[12.6695 0 -1.2353 0 0;
0 1 0 0 0;
0 -1 .36 0 0;
0 0 1 0 0;
0 0 0 1 -.03475];

A=[12.353 0 0 -0.9186 0;
0 .95 0 0 0;
.36 0 0 -0.64 0;
1 0 0 0 1;
0 0 0 1 0];

% Q*A*Z==AA, Q*B*Z==BB
[AA,BB,Q,Z,V,W]=qz(A,B);

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享1 收藏1 回帖

关键词：MATLAB程序 MATLAB atlab matla ABCs matlab 绝对值 程序

问题已经自行解决！用了Sim C. 的程序

呵呵，用几个简单命令也可以。先求特征根和伴随矩阵，然后用SORT命令按特征根的绝对值排序，同时将伴随特征向量也从新排序，就可以了。当然，用Sims的程序也很好。

lustboy 发表于 2012-8-15 20:39
呵呵，用几个简单命令也可以。先求特征根和伴随矩阵，然后用SORT命令按特征根的绝对值排序，同时将伴随特征 ...

不知道你有没有亲自实践过？如果矩阵不能对角化，似乎不能做啊

ahnulxy 发表于 2012-8-16 10:35
不知道你有没有亲自实践过？如果矩阵不能对角化，似乎不能做啊

你是对的，当矩阵正定时比较好用，如果非正定，则比较复杂。
一般情况下，特别是对于纯线性系统，我用上面说的方法，特殊的时候我是用Sargent的程序。

lustboy 发表于 2012-8-16 17:34
你是对的，当矩阵正定时比较好用，如果非正定，则比较复杂。
一般情况下，特别是对于纯线性系统，我用上 ...

谢谢！还没用过Sargent的程序，能告诉在哪能下载到嘛

ahnulxy 发表于 2012-8-16 19:31
谢谢！还没用过Sargent的程序，能告诉在哪能下载到嘛

Sargent个人主页上就有。原来他在Stanford的时候有FTP，里面有很多东西，现在没有了。
ZhaTao那里也有很多程序，他跟Sims合作很多，跟Sargent也有合作。

还有，我也用这个：
% Written by: Yunkai Zhou (周云开), yzhou@smu.edu
% Assistant Professor, Department of Mathematics, Southern Methodist University
% 2005,06,25
%
% this script is adapted from Gerard Sleijpen's jdqr.m
%
function [Q, S, indx]=sortschur_old(Q, S, sigma, kk)

  n = size(S,1);
  
  if n<2,
    indx=1; return,
  elseif nargin<4, kk=n-1; end
  
  if nargin<=4,
    kk=min(kk,n-1);
  else
    sigma=sigma(1,:);
  end

  % transform real schur form to complex schur form if necessary
  if  norm(tril(S,-1),1)>0, [Q,S]=rsf2csf(Q, S); end

  %------ find order eigenvalues ---------------
  indx = SortEig(diag(S),sigma);

  %------ reorder schur form -------------------
  [Q, S] = SwapSchur(Q, S, indx(1:kk));

  return
  
%----------------------------------------------------------------------
function I=SortEig(t,sigma);
%
% I=SortEig(T,SIGMA) sorts T according to SIGMA.
%
% T is a vector of scalars,
% SIGMA is a string or a vector of scalars.
% I is a permutation of (1:LENGTH(T))' such that:
%   if SIGMA is a vector of scalars then
%   for K=1,2,...,LENGTH(T) with KK = MIN(K,SIZE(SIGMA,1))
%      ABS( T(I(K))-SIGMA(KK) ) <= ABS( T(I(J))-SIGMA(KK) )
%      SIGMA(kk)=INF: ABS( T(I(K)) ) >= ABS( T(I(J)) )
%         for all J >= K

if ischar(sigma)
  switch sigma
   case 'LM'
    [s,I]=sort(-abs(t));
   case 'SM'
    [s,I]=sort(abs(t));
   case 'LR'
    [s,I]=sort(-real(t));
   case 'SR'
    [s,I]=sort(real(t));
   case 'BE'
    [s,I]=sort(real(t)); I=twistdim(I,1);
   case 'LI'
    [s,I]=sort(-imag(t));
   case 'SI'
    [s,I]=sort(imag(t));   
  end
else

  [s,I]=sort(abs(t-sigma(1,1)));
  ll=min(size(sigma,1),size(t,1)-1);
  for j=2:ll
    if sigma(j,1)==inf
      [s,J]=sort(abs(t(I(j:end)))); J=flipdim(J,1);
    else
      [s,J]=sort(abs(t(I(j:end))-sigma(j,1)));
    end
    I=[I(1:j-1);I(J+j-1)];
  end

end

return

%----------------------------------------------------------------------
function t=twistdim(t,k)

  d=size(t,k); J=1:d; J0=zeros(1,2*d);
  J0(1,2*J)=J; J0(1,2*J-1)=flipdim(J,2); I=J0(1,J);
  if k==1, t=t(I,:); else, t=t(:,I); end

return

%----------------------------------------------------------------------
function [Q,S]=SwapSchur(Q,S,I)
% [Q,S]=SwapSchur(QQ,SS,P)
%    QQ and SS are square matrices of size K by K
%    P is the first part of a permutation of (1:K)'.
%
%    If    M = QQ*SS*QQ'  and  QQ'*QQ = EYE(K), SS upper triangular
%    then  M*Q = Q*S      with   Q'*Q = EYE(K),  S upper triangular
%    and   D(1:LENGTH(P))=DD(P) where D=diag(S), DD=diag(SS)
%
%    Computation uses Givens rotations. (similar to G.V. 1996, p. 366)

%
% note that using SwapSchur for symmetric problems is not necessary,
% symmetric problems only require direct swap by assignment
%
  
  kk=min(length(I),size(S,1)-1);
  j=1; while (j<=kk & j==I(j)), j=j+1; end;
  while j<=kk
    i=I(j);
    for k=i-1:-1:j
      q = [S(k,k)-S(k+1,k+1),S(k,k+1)];
      if q(1) ~= 0
        q = q/norm(q);
        G = [[q(2);-q(1)],q'];
        J = [k,k+1];
        Q(:,J) = Q(:,J)*G;
        S(:,J) = S(:,J)*G;  
        S(J,:) = G'*S(J,:);
      end
      S(k+1,k) = 0;
    end
    I=I+(I<i);
    j=j+1; while (j<=kk & j==I(j)), j=j+1; end
  end

return

lustboy 发表于 2012-8-16 20:21
Sargent个人主页上就有。原来他在Stanford的时候有FTP，里面有很多东西，现在没有了。
ZhaTao那里也有很 ...

谢谢

学习一下。我也正打算看这本书呢。