octave_packages/m/special-matrix/hilb.m

   1 ## Copyright (C) 1993-2012 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 3 of the License, or (at
   8 ## your option) 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, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {} hilb (@var{n})
  21 ## Return the Hilbert matrix of order @var{n}.  The @math{i,j} element
  22 ## of a Hilbert matrix is defined as
  23 ## @tex
  24 ## $$
  25 ## H(i, j) = {1 \over (i + j - 1)}
  26 ## $$
  27 ## @end tex
  28 ## @ifnottex
  29 ##
  30 ## @example
  31 ## H(i, j) = 1 / (i + j - 1)
  32 ## @end example
  33 ##
  34 ## @end ifnottex
  35 ##
  36 ## Hilbert matrices are close to being singular which make them difficult to
  37 ## invert with numerical routines.
  38 ## Comparing the condition number of a random matrix 5x5 matrix with that of
  39 ## a Hilbert matrix of order 5 reveals just how difficult the problem is.
  40 ##
  41 ## @example
  42 ## @group
  43 ## cond (rand (5))
  44 ##    @result{} 14.392
  45 ## cond (hilb (5))
  46 ##    @result{} 4.7661e+05
  47 ## @end group
  48 ## @end example
  49 ##
  50 ## @seealso{invhilb}
  51 ## @end deftypefn
  52
  53 ## Author: jwe
  54
  55 function retval = hilb (n)
  56
  57   if (nargin != 1)
  58     print_usage ();
  59   elseif (! isscalar (n))
  60     error ("hilb: N must be a scalar integer");
  61   endif
  62
  63   retval = zeros (n);
  64   tmp = 1:n;
  65   for i = 1:n
  66     retval(i, :) = 1.0 ./ tmp;
  67     tmp++;
  68   endfor
  69
  70 endfunction
  71
  72
  73 %!assert (hilb (2), [1, 1/2; 1/2, 1/3])
  74 %!assert (hilb (3), [1, 1/2, 1/3; 1/2, 1/3, 1/4; 1/3, 1/4, 1/5])
  75
  76 %!error hilb ()
  77 %!error hilb (1, 2)
  78 %!error <N must be a scalar integer> hilb (ones(2))
  79