Mercurial > hg > octave-lyh
annotate scripts/special-matrix/invhilb.m @ 11188:4cb1522e4d0f
Use function handle as input to cellfun,
rather than quoted function name or anonymous function wrapper.
| author | Rik <octave@nomad.inbox5.com> |
|---|---|
| date | Wed, 03 Nov 2010 17:20:56 -0700 |
| parents | a4f482e66b65 |
| children | fd0a3ac60b0e |
| rev | line source |
|---|---|
| 7017 | 1 ## Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2002, |
| 8920 | 2 ## 2004, 2005, 2006, 2007, 2008, 2009 Dirk Laurie |
| 2313 | 3 ## |
| 4 ## This file is part of Octave. | |
| 5 ## | |
| 6 ## Octave is free software; you can redistribute it and/or modify it | |
| 7 ## under the terms of the GNU General Public License as published by | |
| 7016 | 8 ## the Free Software Foundation; either version 3 of the License, or (at |
| 9 ## your option) any later version. | |
| 2313 | 10 ## |
| 11 ## Octave is distributed in the hope that it will be useful, but | |
| 12 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
| 13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
| 14 ## General Public License for more details. | |
| 15 ## | |
| 16 ## You should have received a copy of the GNU General Public License | |
| 7016 | 17 ## along with Octave; see the file COPYING. If not, see |
| 18 ## <http://www.gnu.org/licenses/>. | |
| 245 | 19 |
| 3369 | 20 ## -*- texinfo -*- |
| 21 ## @deftypefn {Function File} {} invhilb (@var{n}) | |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
22 ## Return the inverse of the Hilbert matrix of order @var{n}. This can be |
| 5435 | 23 ## computed exactly using |
| 3889 | 24 ## @tex |
| 25 ## $$\eqalign{ | |
| 26 ## A_{ij} &= -1^{i+j} (i+j-1) | |
| 27 ## \left( \matrix{n+i-1 \cr n-j } \right) | |
| 28 ## \left( \matrix{n+j-1 \cr n-i } \right) | |
| 29 ## \left( \matrix{i+j-2 \cr i-2 } \right)^2 \cr | |
| 30 ## &= { p(i)p(j) \over (i+j-1) } | |
| 31 ## }$$ | |
| 32 ## where | |
| 33 ## $$ | |
| 34 ## p(k) = -1^k \left( \matrix{ k+n-1 \cr k-1 } \right) | |
| 35 ## \left( \matrix{ n \cr k } \right) | |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
36 ## $$ |
| 3889 | 37 ## @end tex |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7411
diff
changeset
|
38 ## @ifnottex |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
39 ## |
| 3889 | 40 ## @example |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
41 ## @group |
| 3889 | 42 ## |
| 43 ## (i+j) /n+i-1\ /n+j-1\ /i+j-2\ 2 | |
| 44 ## A(i,j) = -1 (i+j-1)( )( ) ( ) | |
| 45 ## \ n-j / \ n-i / \ i-2 / | |
| 46 ## | |
| 47 ## = p(i) p(j) / (i+j-1) | |
| 48 ## | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
49 ## @end group |
| 3889 | 50 ## @end example |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
51 ## |
|
10846
a4f482e66b65
Grammarcheck more of the documentation.
Rik <octave@nomad.inbox5.com>
parents:
10821
diff
changeset
|
52 ## @noindent |
| 3889 | 53 ## where |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
54 ## |
| 3889 | 55 ## @example |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
56 ## @group |
| 3889 | 57 ## k /k+n-1\ /n\ |
| 58 ## p(k) = -1 ( ) ( ) | |
| 59 ## \ k-1 / \k/ | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
60 ## @end group |
| 3889 | 61 ## @end example |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
62 ## |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7411
diff
changeset
|
63 ## @end ifnottex |
| 3889 | 64 ## The validity of this formula can easily be checked by expanding |
| 65 ## the binomial coefficients in both formulas as factorials. It can | |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
66 ## be derived more directly via the theory of Cauchy matrices. |
|
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
67 ## See J. W. Demmel, @cite{Applied Numerical Linear Algebra}, p. 92. |
| 3889 | 68 ## |
| 69 ## Compare this with the numerical calculation of @code{inverse (hilb (n))}, | |
| 3369 | 70 ## which suffers from the ill-conditioning of the Hilbert matrix, and the |
| 71 ## finite precision of your computer's floating point arithmetic. | |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
72 ## @seealso{hilb, hankel, vander, sylvester_matrix, toeplitz} |
| 3369 | 73 ## @end deftypefn |
| 4 | 74 |
| 5132 | 75 ## Author: Dirk Laurie <dlaurie@na-net.ornl.gov> |
| 2314 | 76 |
| 2311 | 77 function retval = invhilb (n) |
| 4 | 78 |
| 79 if (nargin != 1) | |
| 6046 | 80 print_usage (); |
| 4 | 81 endif |
| 82 | |
| 83 nmax = length (n); | |
| 84 if (nmax == 1) | |
| 3889 | 85 |
| 86 ## The point about the second formula above is that when vectorized, | |
| 87 ## p(k) is evaluated for k=1:n which involves O(n) calls to bincoeff | |
| 88 ## instead of O(n^2). | |
| 89 ## | |
| 90 ## We evaluate the expression as (-1)^(i+j)*(p(i)*p(j))/(i+j-1) except | |
| 91 ## when p(i)*p(j) would overflow. In cases where p(i)*p(j) is an exact | |
| 92 ## machine number, the result is also exact. Otherwise we calculate | |
| 93 ## (-1)^(i+j)*p(i)*(p(j)/(i+j-1)). | |
| 94 ## | |
| 4031 | 95 ## The Octave bincoeff routine uses transcendental functions (gammaln |
| 3889 | 96 ## and exp) rather than multiplications, for the sake of speed. |
| 97 ## However, it rounds the answer to the nearest integer, which | |
| 98 ## justifies the claim about exactness made above. | |
| 99 | |
| 100 retval = zeros (n); | |
| 101 k = [1:n]; | |
| 102 p = k .* bincoeff (k+n-1, k-1) .* bincoeff (n, k); | |
| 103 p(2:2:n) = -p(2:2:n); | |
| 104 if (n < 203) | |
| 105 for l = 1:n | |
| 106 retval(l,:) = (p(l) * p) ./ [l:l+n-1]; | |
| 4 | 107 endfor |
| 3889 | 108 else |
| 109 for l = 1:n | |
| 110 retval(l,:) = p(l) * (p ./ [l:l+n-1]); | |
| 111 endfor | |
| 112 endif | |
| 4 | 113 else |
| 3889 | 114 error ("invhilb: expecting scalar argument, found something else"); |
| 4 | 115 endif |
| 116 | |
| 117 endfunction | |
| 7411 | 118 |
| 119 %!test | |
| 120 %! result4 = [16, -120, 240, -140; | |
| 121 %! -120, 1200, -2700, 1680; | |
| 122 %! 240, -2700, 6480, -4200; | |
| 123 %! -140, 1680, -4200, 2800]; | |
| 124 %! | |
| 125 %! assert((invhilb (1) == 1 && invhilb (2) == [4, -6; -6, 12] | |
| 126 %! && invhilb (4) == result4 | |
| 127 %! && abs (invhilb (7) * hilb (7) - eye (7)) < sqrt (eps))); | |
| 128 | |
| 129 %!error invhilb ([1, 2]); | |
| 130 | |
| 131 %!error invhilb (); | |
| 132 | |
| 133 %!error invhilb (1, 2); | |
| 134 |