home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / r / rlab / CTB / lyap2 < prev    next >
Encoding:
Text File  |  1995-11-15  |  2.5 KB  |  89 lines
  1. //-----------------------------------------------------------------------
  2. //
  3. // lyap2
  4. //
  5. // SyntaX: x=lyap2(A,B,C)
  6. //
  7. // This routine solves the lyapunov equation using an eigenvalue
  8. // decomposition (spectral decomposition). The solution to the
  9. // equation is returned in a matrix X (X=lyap2(A,B,C) ).
  10. //
  11. // If the routine is called with, X=lyap2(A,B) then it solves the
  12. // particular form of the Lyapunov equation:
  13. //
  14. //       A*X+X*A' = -B
  15. //
  16. // Calling the routine with X=lyap(A,B,C) solves the general form
  17. // of the lyapunov equation (sylvester equation):
  18. //
  19. //       A*X+X*B = -C
  20. //
  21. // ==================================================================
  22. // Note: This routine assumes that there are no repeated eigenvalues.
  23. //       If repeated eigenvalues exist, then the solution may be
  24. //       unreliable.
  25. // ==================================================================
  26. //
  27. // Ref: Junkins, J., and Kim, Y., "Introduction to Dynamics and Control
  28. //      of Flexible Structures," AIAA Inc, Washington D.C., 1993,
  29. //      Chapter 2.
  30. //
  31. // Copyright (C), by Jeffrey B. Layton, 1994
  32. // Version JBL 940919
  33. //-----------------------------------------------------------------------
  34.  
  35. lyap2 = function(A,B,C)
  36. {
  37.    local(nargs,PhiA,LamA,PhiB,LamB,E,X,PsiA,PsiB,T,DA,DB);
  38.  
  39. // Count number of input arguments
  40.    nargs=0;
  41.    if (exist(A)) {nargs=nargs+1;}
  42.    if (exist(B)) {nargs=nargs+1;}
  43.    if (exist(C)) {nargs=nargs+1;}
  44.  
  45.    if (nargs != 3) {
  46.        error("LYAP2: Wrong number of input arguments.");
  47.    }
  48.  
  49. // Dimension Check
  50.    if (A.nr != A.nc) { error("lyap2: Dimensions of A and B do not agree"); }
  51.    if (B.nr != B.nc) { error("lyap2: Dimensions of A and B do not agree"); }
  52.    if (A.nr == 0) {
  53.        X = [];
  54.        return X;
  55.    }
  56.  
  57. // Compute eigenvalues and the left and right eigenvectors of A
  58.    E=eig(A);
  59.    LamA=diag(E.val);
  60.    DA=LamA;
  61.    PhiA=E.vec;
  62.    PsiA=eig(A').vec;
  63.  
  64. // Switch on Lyapunov solution or Sylvester solution.
  65.    if ( exist (C)) {
  66.        if (C.nr != A.nr) { error("lyap2: Dimensions of C do not agree"); }
  67.        if (C.nc != B.nr) { error("lyap2: Dimensions of C do not agree"); }
  68.  
  69. // Computes eigenvalues and left and right eigenvectors of B
  70.        E=eig(B);
  71.        LamB=diag(E.val);
  72.        DB=LamB;
  73.        LamB=ones(A.nr,B.nr)*LamB;
  74.        PhiB=E.vec;
  75.        PsiB=eig(B').vec;
  76.        LamA=LamA*ones(A.nr,B.nr);
  77.  
  78.        T=-(PsiA.'*C*PhiB) ./ (LamA+LamB);
  79.        X=PhiA*T*PhiB';
  80.    else
  81.        LamA=LamA*ones(A.nr,A.nr);
  82.        T=-(PsiA.'*B*PsiA) ./ (LamA+LamA.');
  83.        X=PhiA*T*PhiA';
  84.    }
  85.    X=real(X);
  86.  
  87.    return X;
  88. };
  89.