home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / hnfsnric.m < prev    next >
Text File  |  1999-04-29  |  2KB  |  73 lines
  1. # Copyright (C) 1996,1998 A. Scottedward Hodel
  2. #
  3. # This file is part of Octave.
  4. #
  5. # Octave is free software; you can redistribute it and/or modify it
  6. # under the terms of the GNU General Public License as published by the
  7. # Free Software Foundation; either version 2, or (at your option) any
  8. # later version.
  9. #
  10. # Octave is distributed in the hope that it will be useful, but WITHOUT
  11. # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12. # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13. # for more details.
  14. #
  15. # You should have received a copy of the GNU General Public License
  16. # along with Octave; see the file COPYING.  If not, write to the Free
  17. # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  
  19. function [Xinf,x_ha_err] = hnfsnric(A,BB,C1,d1dot,R,ptol)
  20. #
  21. # forms
  22. #        xx = ([BB; -C1'*d1dot]/R) * [d1dot'*C1 BB'];
  23. #        Ha = [A 0*A; -C1'*C1 -A'] - xx;
  24. # and solves associated Riccati equation
  25. # returns error code
  26. #  x_ha_err:
  27. #    0: successful
  28. #    1: Xinf has imaginary eigenvalues
  29. #    2: Hx not Hamiltonian
  30. #    3: Xinf has inf. eigenvalues (numerical overflow)
  31. #    4: Xinf not symmetric
  32. #    5: Xinf not positive definite
  33. #    6: R is singular
  34.  
  35. x_ha_err = 0;        # assume success
  36. Xinf = [];             # default return value
  37. n = is_sqr(A);
  38. nw = is_sqr(R);
  39. if(rank(R) != nw)    x_ha_err = 6;
  40. else                 # build hamiltonian Ha for X_inf
  41.   xx = ([BB; -C1'*d1dot]/R) * [d1dot'*C1, BB'];
  42.   Ha = [A, 0*A; -C1'*C1, -A'] - xx;
  43.   x_ha_err = 0;
  44.   [d, Ha] = balance(Ha);
  45.   [u, s] = schur(Ha, "A");
  46.   rev = real(eig(s));
  47.  
  48.   if (any(abs(rev) <= ptol))    # eigenvalues near the imaginary axis
  49.     x_ha_err = 1;
  50.   elseif (sum(rev > 0) != sum(rev < 0))
  51.     # unequal number of positive and negative eigenvalues
  52.     x_ha_err = 2;
  53.   else
  54.     # compute positive Riccati equation solution
  55.     u = d * u;
  56.     Xinf = u(n+1:2*n,1:n) / u(1:n,1:n);
  57.     if (!all(all(finite(Xinf))))
  58.       x_ha_err = 3;
  59.     elseif (norm(Xinf-Xinf') >= 10*ptol)
  60.       # solution not symmetric
  61.       x_ha_err = 4;
  62.     else
  63.       # positive semidefinite?
  64.       # force symmetry (faster, avoids some convergence problems)
  65.       Xinf = (Xinf + Xinf')/2;
  66.       rev = eig(Xinf);
  67.       if (any(rev <= -ptol))
  68.         x_ha_err = 5;
  69.       endif
  70.     endif
  71.   endif
  72. endif
  73.