home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / hnfsnric.m < prev    next >
Text File  |  1999-12-24  |  3KB  |  76 lines
  1. ## Copyright (C) 1996,1998 Auburn University.  All Rights Reserved
  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, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
  18.  
  19. ## usage: [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. function [Xinf, x_ha_err] = hnfsnric (A, BB, C1, d1dot, R, ptol)
  36.  
  37.   x_ha_err = 0;        # assume success
  38.   Xinf = [];             # default return value
  39.   n = is_sqr(A);
  40.   nw = is_sqr(R);
  41.   if(rank(R) != nw)    x_ha_err = 6;
  42.   else                 # build hamiltonian Ha for X_inf
  43.     xx = ([BB; -C1'*d1dot]/R) * [d1dot'*C1, BB'];
  44.     Ha = [A, 0*A; -C1'*C1, -A'] - xx;
  45.     x_ha_err = 0;
  46.     [d, Ha] = balance(Ha);
  47.     [u, s] = schur(Ha, "A");
  48.     rev = real(eig(s));
  49.  
  50.     if (any(abs(rev) <= ptol))    # eigenvalues near the imaginary axis
  51.       x_ha_err = 1;
  52.     elseif (sum(rev > 0) != sum(rev < 0))
  53.       ## unequal number of positive and negative eigenvalues
  54.       x_ha_err = 2;
  55.     else
  56.       ## compute positive Riccati equation solution
  57.       u = d * u;
  58.       Xinf = u(n+1:2*n,1:n) / u(1:n,1:n);
  59.       if (!all(all(finite(Xinf))))
  60.     x_ha_err = 3;
  61.       elseif (norm(Xinf-Xinf') >= 10*ptol)
  62.     ## solution not symmetric
  63.     x_ha_err = 4;
  64.       else
  65.     ## positive semidefinite?
  66.     ## force symmetry (faster, avoids some convergence problems)
  67.     Xinf = (Xinf + Xinf')/2;
  68.     rev = eig(Xinf);
  69.     if (any(rev <= -ptol))
  70.       x_ha_err = 5;
  71.     endif
  72.       endif
  73.     endif
  74.   endif
  75. endfunction
  76.