home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / orthog_r < prev    next >
Encoding:
Text File  |  1994-12-20  |  3.3 KB  |  105 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Orthogonal and nearly orthogonal matrices.
  4.  
  5. // Syntax:      Q = orthog ( N , K )
  6.  
  7. // Description:
  8.  
  9. //      Orthog produces the K'th type of matrix of order N.
  10.  
  11. //      K > 0 for exactly orthogonal matrices, K < 0 for diagonal
  12. //      scalings of orthogonal matrices.
  13.  
  14. //      Available types: (K = 1 is the default)
  15. //       K = 1:  Q[i;j] = SQRT(2/(n+1)) * SIN( i*j*PI/(n+1) )
  16. //               Symmetric eigenvector matrix for second difference matrix.
  17. //       K = 2:  Q[i;j] = 2/SQRT(2*n+1)) * SIN( 2*i*j*PI/(2*n+1) )
  18. //               Symmetric.
  19. //       K = 3:  Q[r;s] = EXP(2*PI*i*(r-1)*(s-1)/n) / SQRT(n)  (i=SQRT(-1))
  20. //               Unitary, the Fourier matrix.  Q^4 is the identity.
  21. //               This is essentially the same matrix as FFT(EYE(N))/SQRT(N)!
  22. //       K = 4:  Helmert matrix: a permutation of a lower Hessenberg matrix,
  23. //               whose first row is ONES(1:N)/SQRT(N).
  24. //       K = 5:  Q[i;j] = SIN( 2*PI*(i-1)*(j-1)/n ) + COS( 2*PI*(i-1)*(j-1)/n ).
  25. //               Symmetric matrix arising in the Hartley transform.
  26. //       K = -1: Q[i;j] = COS( (i-1)*(j-1)*PI/(n-1) )
  27. //               Chebyshev Vandermonde-like matrix, based on extrema of T(n-1).
  28. //       K = -2: Q[i;j] = COS( (i-1)*(j-1/2)*PI/n) )
  29. //               Chebyshev Vandermonde-like matrix, based on zeros of T(n).
  30.  
  31. //      References:
  32. //       N.J. Higham and D.J. Higham, Large growth factors in Gaussian
  33. //            elimination with pivoting, SIAM J. Matrix Analysis and  Appl.,
  34. //            10 (1989), pp. 155-164.
  35. //       P. Morton, On the eigenvectors of Schur's matrix, J. Number Theory,
  36. //            12 (1980), pp. 122-127. (Re. ORTHOG(N, 3))
  37. //       H.O. Lancaster, The Helmert Matrices, Amer. Math. Monthly, 72 (1965),
  38. //            pp. 4-12.
  39. //       D. Bini and P. Favati, On a matrix algebra related to the discrete
  40. //            Hartley transform, SIAM J. Matrix Anal. Appl., 14 (1993),
  41. //            pp. 500-507.
  42.  
  43. //    This file is a translation of orthog.m from version 2.0 of
  44. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  45. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  46.  
  47. //-------------------------------------------------------------------//
  48.  
  49. orthog = function ( n , k )
  50. {
  51.   local (n, k)
  52.   global (pi)
  53.  
  54.   if (!exist (k)) { k = 1; }
  55.  
  56.   if (k == 1)
  57.   {
  58.     // E'vectors second difference matrix
  59.     m = (1:n)'*(1:n) * (pi/(n+1));
  60.     Q = sin(m) * sqrt(2/(n+1));
  61.  
  62.   else if (k == 2) {
  63.  
  64.     m = (1:n)'*(1:n) * (2*pi/(2*n+1));
  65.     Q = sin(m) * (2/sqrt(2*n+1));
  66.  
  67.   else if (k == 3) {
  68.     //  Vandermonde based on roots of unity
  69.  
  70.     m = 0:n-1;
  71.     Q = exp(m'*m*2*pi*sqrt(-1)/n) / sqrt(n);
  72.  
  73.   else if (k == 4) {
  74.     //  Helmert matrix
  75.     Q = tril(ones(n,n));
  76.     Q[1;2:n] = ones(1,n-1);
  77.     for (i in 2:n)
  78.     {
  79.       Q[i;i] = -(i-1);
  80.     }
  81.     Q = diag( sqrt( [n, 1:n-1] .* [1:n] ) ) \ Q;
  82.  
  83.   else if (k == 5) {
  84.     //  Hartley matrix
  85.     m = (0:n-1)'*(0:n-1) * (2*pi/n);
  86.     Q = (cos(m) + sin(m))/sqrt(n);
  87.  
  88.   else if (k == -1) {
  89.     //  extrema of T(n-1)
  90.     m = (0:n-1)'*(0:n-1) * (pi/(n-1));
  91.     Q = cos(m);
  92.  
  93.   else if (k == -2) {
  94.     // zeros of T(n)
  95.     m = (0:n-1)'*(.5:n-.5) * (pi/n);
  96.     Q = cos(m);
  97.  
  98.   else
  99.     error ("Illegal value for second parameter.");
  100.  
  101.   } } } } } } }
  102.  
  103.   return Q;
  104. };
  105.