home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / r / rlab / CTB / eigsen2 < prev    next >
Encoding:
Text File  |  1995-11-15  |  3.3 KB  |  111 lines
  1. //-------------------------------------------------------------------------------
  2. //
  3. // eigsen
  4. //
  5. // Syntax: D=eigsen(A,B,DADP,DBDP)
  6. //
  7. // This routine computes the gradient of the eigenvalues, right eigenvectors,
  8. // and left eigenvectors of a matrix with respect to a scalar parameter.It
  9. // considers a general eigenvalue problem of the form:
  10. //
  11. //    Ax = (lambda) B x
  12. //
  13. // with the normalizations:
  14. //
  15. //    phi(i)'*B*phi(i) = 1
  16. //    psi(i)'*A*phi(j) = krnecker delta(i,j)
  17. //
  18. // where phi is right eigenvector matrix, psi is the let eigenvector matrix,
  19. // and where phi(i) is the ith right eigenvector and psi(i) is the ith left
  20. // eigenvector, krnecker delta is the Kronecker Delta function.
  21. //
  22. // If the routine is called as,
  23. //
  24. //    D=eigsen2(A,B,DADP,DBDP)
  25. //
  26. // where DADP and DBDP are the gradients of A and B respectively, with respect
  27. // to a scalar p, it returns the gradient of the eigenvalues (D.de), the
  28. // gradient of the right eigenvectors (D.dPhi), and the gradient of the
  29. // left eigenvectors (D.dPsi).
  30. //
  31. // The results are returned in a list:
  32. //
  33. //     D.de = Gradient of the Eigenvalues
  34. //   D.dPhi = Gradient of the Right Eigenvectors
  35. //   D.dPsi = Gradient of the Left Eigenvectors
  36. //
  37. // Note: The routine makes an effort to check if there are no repeated
  38. //       eigenvalues. However, only a warning is printed and the routine
  39. //       continues. It is up to the user to make sure there are no repeated
  40. //       eigenvalues.
  41. //
  42. // Refs: (1) Lim, K.B., Junkins, J.L., and Wang, B.P.
  43. //           "Re-examination of Eigenvector Derivatives"
  44. //           Journal of Guidance, Control, and Dynamics, Vol. 10, No. 6, 1987,
  45. //           pp. 581-587.
  46. //
  47. //       (2) Junkins and Kim, Dyn. & Ctrl of Structures, Introduction to Dynamics
  48. //           and Control of Flexible Structures, AIAA, 1993.
  49. //
  50. // Copyright (C), by Jeffrey B. Layton, 1994
  51. // Version JBL 940918
  52. //-------------------------------------------------------------------------------
  53.  
  54. rfile eignm
  55.  
  56. eigsen2 = function(A,B,DADP,DBDP)
  57. {
  58.    local(n,a,b,xi,yi,xk,yk,de,dPhi,dPsi,ssum,i,k,Adum,Lambda,nargs)
  59.  
  60. // Count number of input arguments
  61.    nargs=0;
  62.    if (exist(A)) {nargs=nargs+1;}
  63.    if (exist(B)) {nargs=nargs+1;}
  64.    if (exist(DADP)) {nargs=nargs+1;}
  65.    if (exist(DBDP)) {nargs=nargs+1;}
  66.  
  67.    if (nargs < 4) {
  68.        error("EIGSEN2: Wrong number of input arguments.");
  69.    }
  70.  
  71. // Compute the spectral decomposition using the correct normalization
  72. // by using the routine "eign"
  73.    Adum=eignm(A,B);
  74.    Lambda=Adum.Lambda;
  75.  
  76. // Initialization
  77.    n=max(size(Lambda));
  78.  
  79. // Initialize Results
  80.    a=zeros(n,1);
  81.    de=zeros(n,1);
  82.    dPhi=zeros(n,n);
  83.    dPsi=zeros(n,n);
  84.  
  85. // Begin Computations
  86.    for (i in 1:n) {
  87.         xi=Adum.Phi[;i];
  88.         yi=Adum.Psi[;i];
  89.         de[i]=yi.'*(DADP-Lambda[i]*DBDP)*xi;
  90.         ssum=0;
  91.         for (k in 1:n) {
  92.              if (k != i) {
  93.                  xk=Adum.Phi[;k];
  94.                  yk=Adum.Psi[;k];
  95.                  a[k]=yk.'*(DADP-Lambda[i]*DBDP)*xi/(Lambda[i]-Lambda[k]);
  96.                  ssum=ssum+a[k]*xk.'*(B+B')*xi;
  97.                  b[k]=yi.'*(DADP-Lambda[i]*DBDP)*xk/(Lambda[i]-Lambda[k]);
  98.              }
  99.         }
  100. // Add in term when k=i
  101.         a[i]=-(xi.'*DBDP*xi+ssum)/2;
  102.         b[i]=-yi.'*DBDP*xi-a[i];
  103.         for (k in 1:n) {
  104.              dPhi[;i]=dPhi[;i]+a[k]*Adum.Phi[;k];
  105.              dPsi[;i]=dPsi[;i]+b[k]*Adum.Psi[;k];
  106.         }
  107.    }
  108.  
  109.    return << de=de; dPhi=dPhi; dPsi=dPsi >>
  110. };
  111.