home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts / special-matrix / invhilb.m < prev    next >
Text File  |  1998-11-10  |  2KB  |  68 lines
  1. ## Copyright (C) 1996, 1997 John W. Eaton
  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
  7. ## the Free Software Foundation; either version 2, or (at your option)
  8. ## any later version.
  9. ##
  10. ## Octave is distributed in the hope that it will be useful, but
  11. ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12. ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13. ## General Public License 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
  18. ## 02111-1307, USA.
  19.  
  20. ## usage: invhilb (n)
  21. ##
  22. ## Return the inverse of a Hilbert matrix of order n.  This is slow but
  23. ## exact.  Compare with inv (hilb (n)).
  24. ##
  25. ## See also: hankel, vander, sylvester_matrix, hilb, toeplitz
  26.  
  27. ## Author: jwe
  28.  
  29. function retval = invhilb (n)
  30.  
  31.   if (nargin != 1)
  32.     usage ("invhilb (n)");
  33.   endif
  34.  
  35.   nmax = length (n);
  36.   if (nmax == 1)
  37.     retval = zeros (n);
  38.     for l = 1:n
  39.       for k = l:n
  40.         tmp = 1;
  41.         for i = 1:n
  42.           tmp = tmp * (i + k - 1);
  43.         endfor
  44.         for i = 1:n
  45.           if (i != k)
  46.             tmp = tmp * (l + i - 1);
  47.           endif
  48.         endfor
  49.         for i = 1:n
  50.           if (i != l)
  51.             tmp = tmp / (i - l);
  52.           endif
  53.         endfor
  54.         for i = 1:n
  55.           if (i != k)
  56.             tmp = tmp / (i - k);
  57.           endif
  58.         endfor
  59.         retval (k, l) = tmp;
  60.         retval (l, k) = tmp;
  61.       endfor
  62.     endfor
  63.   else
  64.     error ("hilb: expecting scalar argument, found something else");
  65.   endif
  66.  
  67. endfunction
  68.