home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / ipjfact_r < prev    next >
Encoding:
Text File  |  1994-12-20  |  1.4 KB  |  68 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    A Hankel matrix with factorial elements.
  4.  
  5. // Syntax:      A = ipjfact ( N, K )
  6.  
  7. // Description:
  8.  
  9. //      A.A is the matrix with
  10.  
  11. //                A[i;j] = (i+j)!    (K = 0, default)
  12. //                A[i;j] = 1/(i+j)!  (K = 1)
  13.  
  14. //      Both are Hankel matrices.
  15. //      The determinant and inverse are known explicitly. 
  16.  
  17. //      A list is returned with elements `A' and `detA'
  18.  
  19. //      Suggested by P. R. Graves-Morris.
  20.  
  21. //      Reference:
  22. //      M.J.C. Gover, The explicit inverse of factorial Hankel matrices,
  23. //      Dept. of Mathematics, University of Bradford, 1993.
  24.  
  25. //    This file is a translation of ipjfact.m from version 2.0 of
  26. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  27. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  28.  
  29. // Dependencies
  30.    require hankel
  31.  
  32. //-------------------------------------------------------------------//
  33.  
  34. ipjfact = function ( n , k )
  35. {
  36.   local (n, k)
  37.  
  38.   if (!exist (k)) { k = 0; }
  39.  
  40.   c = cumprod(2:n+1);
  41.   d = cumprod(n+1:2*n) * c[n-1];
  42.  
  43.   A = hankel(c, d);
  44.  
  45.   if (k == 1)
  46.   {
  47.     A = ones(n, n)./A;
  48.   }
  49.  
  50.   d = 1;
  51.   if (k == 0)
  52.   {
  53.     for (i in 1:n-1)
  54.     {
  55.       d = d*prod(1:i+1)*prod(1:n-i);
  56.     }
  57.     d = d*prod(1:n+1);
  58.   else
  59.     for (i in 0:n-1)
  60.     {
  61.       d = d*prod(1:i)/prod(1:n+1+i);
  62.     }
  63.     if (mod(n*(n-1)/2,2)) { d = -d; }
  64.   }
  65.  
  66.   return << A = A; detA = d >>;
  67. };
  68.