home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / ps_r < prev    next >
Encoding:
Text File  |  1994-12-26  |  3.0 KB  |  116 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Dot plot of a pseudospectrum.
  4.  
  5. // Syntax:      ps (A, M, TOL, RL) 
  6.  
  7. // Description:
  8.  
  9. //      Plots an approximation to a pseudospectrum of the square
  10. //      matrix A, using M random perturbations of size TOL.
  11.  
  12. //      M defaults to a SIZE(A)-dependent value and TOL to 1E-3.
  13.  
  14. //      RL defines the type of perturbation:
  15. //         RL =  0 (default): absolute complex perturbations of 2-norm TOL.
  16. //         RL =  1:           absolute real perturbations of 2-norm TOL.
  17. //         RL = -1:           componentwise real perturbations of size TOL.
  18.  
  19. //      The eigenvalues of A are plotted as symbols.
  20.  
  21. //      ps (A, M, TOL, RL, MARKER) creates a plot if MARKER is >= to
  22. //      0. If MARKER is < 0, the no plot is created, but the plot data
  23. //      is returned.
  24.  
  25. //      ps (A, 0) plots just the eigenvalues of A.
  26.  
  27. //      For a given TOL, the pseudospectrum of A is the set of
  28. //      pseudo-eigenvalues of A, that is, the set { e : e is an
  29. //      eigenvalue of A+E, for some E with NORM(E) <= TOL }.
  30. //
  31. //      Reference:
  32. //      L.N. Trefethen, Pseudospectra of matrices, in D.F. Griffiths and
  33. //           G.A. Watson, eds, Numerical Analysis 1991, Proceedings of the 14th
  34. //           Dundee Conference, vol. 260, Pitman Research Notes in Mathematics,
  35. //           Longman Scientific and Technical, Essex, UK, 1992, pp. 234-266.
  36.  
  37. //    This file is a translation of ps.m from version 2.0 of
  38. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  39. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  40.  
  41. // Dependencies
  42.    require cpltaxes
  43.  
  44. //-------------------------------------------------------------------//
  45.  
  46. ps = function (A, m, tol, rl, marksize)
  47. {
  48.   local (A, m, tol, rl, marksize)
  49.  
  50.   if (diff (size (A))) { error ("Matrix must be square."); }
  51.   n = max (size (A));
  52.  
  53.   if (!exist (marksize)) { marksize = 0; }
  54.   if (!exist (rl)) { rl = 0; }
  55.   if (!exist (tol)) { tol = 1e-3; }
  56.   if (!exist (m)) { m = max(1, round( 25*exp(-0.047*n) )); }
  57.  
  58.   if (m == 0)
  59.   {
  60.     e = eig (A).val.';
  61.     ax = cpltaxes (e);
  62.     plimits (ax[1], ax[2], ax[3], ax[4]);
  63.     plstyle ("point"); plegend ();
  64.     plaspect (1); plaxis ();
  65.     plot ([ real (e), imag (e) ]);
  66.     return 1;
  67.   }
  68.  
  69.   x = zeros (m*n,1);
  70.   i = sqrt (-1);
  71.  
  72.   for (j in 1:m)
  73.   {
  74.     if (rl == -1)
  75.     {
  76.       // Componentwise.
  77.       rand ("uniform", -1, 1);
  78.       dA = -ones(n) + 2*rand(n,n);   // Uniform random numbers on [-1,1].
  79.       dA = tol * A .* dA;
  80.     else
  81.       rand ("normal", 0, 1);
  82.       if (rl == 0)
  83.       {
  84.     // Complex absolute.
  85.     dA = rand (n,n) + i*rand (n,n);
  86.       else         
  87.     // Real absolute.
  88.     dA = rand (n,n);
  89.       }
  90.       dA = tol/norm (dA)*dA;
  91.     }
  92.     e = eig(A + dA).val.';
  93.     x[(j-1)*n+1:j*n] = e;
  94.   }
  95.  
  96.   if (marksize >= 0)
  97.   {
  98.     ax = cpltaxes(x);
  99.     plimits (ax[1], ax[2], ax[3], ax[4]);
  100.     plstyle ("point"); plegend ();
  101.     plaspect (1); plaxis ();
  102.     plot ([ real(x), imag(x)] );
  103.  
  104.     e = eig(A).val.';
  105.     subplot (0);
  106.     plot ([ real (e), imag (e) ]);
  107.  
  108.     return x;
  109.  
  110.   else
  111.  
  112.     return x;
  113.  
  114.   }
  115. };
  116.