home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / fv_r < prev    next >
Encoding:
Text File  |  1994-12-27  |  3.7 KB  |  126 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Field of values (or numerical range)
  4.  
  5. // Syntax:      fv (A, NK, THMAX) 
  6.  
  7. // Description:
  8.  
  9. //      fv evaluates and plots the field of values of the NK largest
  10. //      leading principal submatrices of A, using THMAX equally spaced
  11. //      angles in the complex plane.  The defaults are NK = 1 and
  12. //      THMAX = 16.  (For a `publication quality' picture, set THMAX
  13. //      higher, say 32.)  The eigenvalues of A are displayed as `x'.
  14.  
  15. //      Alternative usage: fv(A, NK, THMAX, 1) suppresses the plot and
  16. //      returns the field of values plot data and A's eigenvalues in a
  17. //      list. Note that norm(F,"i") approximates the numerical radius,
  18.  
  19. //              max {abs(z): z is in the field of values of A}.
  20.  
  21. //      Theory:
  22.  
  23. //      Field of values fv(A) = set of all Rayleigh quotients. fv(A)
  24. //      is a convex set containing the eigenvalues of A.  When A is
  25. //      normal fv(A) is the convex hull of the eigenvalues of A (but
  26. //      not vice versa).
  27.  
  28. //               z = x'Ax/(x'x),  z' = x'A'x/(x'x)
  29. //               => real(z) = x'Hx/(x'x),   H = (A+A')/2
  30. //       so      min(eig(H)) <= real(z) <= max(eig(H))
  31.  
  32. //      with equality for x = corresponding eigenvectors of H.  For
  33. //      these x, rq(A,x) is on the boundary of fv(A).
  34.  
  35. //      Based on an original routine by A. Ruhe.
  36. //
  37. //       References:
  38. //       R.A. Horn and C.R. Johnson, Topics in Matrix Analysis, Cambridge
  39. //            University Press, 1991, Section 1.5.
  40. //       A.S. Householder, The Theory of Matrices in Numerical Analysis,
  41. //            Blaisdell, New York, 1964, Section 3.3.
  42. //       C.R. Johnson, Numerical determination of the field of values of a
  43. //            general complex matrix, SIAM J. Numer. Anal., 15 (1978),
  44. //            pp. 595-602.
  45.  
  46. //    This file is a translation of fv.m from version 2.0 of
  47. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  48. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  49.  
  50. // Dependencies
  51.    require rq symmpart skewpart cpltaxes
  52.  
  53. //-------------------------------------------------------------------//
  54.  
  55. fv = function (B, nk, thmax, noplot)
  56. {
  57.   local (B, nk, thmax, noplot)
  58.   global (pi)
  59.  
  60.   if (!exist (nk)) { nk = 1; }
  61.   if (!exist (thmax)) { thmax = 16; }
  62.  
  63.   // Because code below uses thmax + 1 angles.
  64.   thmax = thmax - 1;
  65.  
  66.   iu = sqrt(-1);
  67.   n = B.nr; 
  68.   p = B.nc;
  69.   if (n != p) { error("Matrix must be square."); }
  70.   f = [];
  71.   z = zeros(2*thmax+1,1);
  72.   e = eig(B).val;
  73.   
  74.   // Filter out cases where B is Hermitian or skew-Hermitian, for efficiency.
  75.  
  76.   if (norm(skewpart(B)) == 0)
  77.   {
  78.     f = [min(e), max(e)];
  79.   else if (norm(symmpart(B)) == 0) {
  80.     e = imag(e);
  81.     f = [min(e), max(e)];
  82.     e = iu*e; 
  83.     f = iu*f;
  84.   else
  85.     for (m in 1:nk)
  86.     {
  87.       ns = n+1-m;
  88.       A = B[1:ns; 1:ns];
  89.       for (i in 0:thmax)
  90.       {
  91.     th = i/thmax*pi;
  92.     Ath = exp(iu*th)*A;               // Rotate A through angle th.
  93.     H = 0.5*(Ath + Ath');             // Hermitian part of rotated A.
  94.     XD = eig(H);
  95.     k = sort(real(XD.val)).ind;
  96.     z[1+i] = rq(A, XD.vec[;k[1]]);       // RQ's of A corr. to eigenvalues of H
  97.     z[1+i+thmax] = rq(A,XD.vec[;k[ns]]); // with smallest/largest real part.
  98.       }
  99.       f = [f; z];
  100.     }
  101.  
  102.     // Next line ensures boundary is `joined up' (needed for orthogonal matrices).
  103.     
  104.     f = [f; f[1;]];
  105.     
  106.   } } 
  107.   
  108.   if (thmax == 0) {  f = e; }
  109.  
  110.   if (!exist (noplot))
  111.   {
  112.     ax = cpltaxes(f);
  113.     plaspect (1);
  114.     plegend ();
  115.     plimits (ax[1], ax[2], ax[3], ax[4]);
  116.     plstyle ("line-point");
  117.     plaxis ();
  118.     plot([real(f), imag(f)]);    // Plot the field of values
  119.     subplot(0);
  120.     plstyle ("point");
  121.     plot([real(e)', imag(e)']);  // Plot the eigenvalues too.
  122.   }
  123.  
  124.   return << f=f ; e=e >>;
  125. };
  126.