home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / tridiag_r < prev    next >
Encoding:
Text File  |  1994-12-20  |  2.0 KB  |  73 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Tridiagonal matrix
  4.  
  5. // Syntax:    T = tridiag ( X , Y , Z )
  6.  
  7. // Description:
  8.  
  9. //    T is the tridiagonal matrix with subdiagonal X, diagonal Y,
  10. //    and superdiagonal Z.  X and Z must be vectors of dimension one
  11. //    less than Y. 
  12.  
  13. //    Alternatively TRIDIAG(N, C, D, E), where C, D, and E are all
  14. //    scalars, yields the Toeplitz tridiagonal matrix of order N
  15. //    with subdiagonal elements C, diagonal elements D, and
  16. //    superdiagonal elements E.   
  17.  
  18. //    This matrix has eigenvalues (Todd 1977)
  19. //        D + 2*SQRT(C*E)*COS(k*PI/(N+1)), k=1:N.
  20. //    tridiag(N) is the same as tridiag(N,-1,2,-1), which is a
  21. //    symmetric positive definite M-matrix (the negative of the
  22. //    second difference matrix). 
  23.  
  24. //      References:
  25. //         J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
  26. //           Birkhauser, Basel, and Academic Press, New York, 1977, p. 155.
  27. //         D.E. Rutherford, Some continuant determinants arising in physics and
  28. //           chemistry---II, Proc. Royal Soc. Edin., 63, A (1952), pp. 232-241.
  29.  
  30. //    This file is a translation of tridiag.m from version 2.0 of
  31. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  32. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  33.  
  34. //-------------------------------------------------------------------//
  35.  
  36. tridiag = function ( n , x , y , z )
  37. {
  38.   local (n, x, y, z)
  39.  
  40.   if (!exist (x))    // nargin = 1
  41.   {
  42.     x = -1; y = 2; z = -1;
  43.   }
  44.  
  45.   if (!exist (z))    // nargin = 3
  46.   {
  47.     z = y; y = x; x = n;
  48.   }
  49.  
  50.   x = x[:]; y = y[:]; z = z[:];        // Force column vectors.
  51.  
  52.   if (max( [ size(x), size(y), size(z) ] ) == 1)
  53.   {
  54.     x = x*ones(n-1,1);
  55.     z = z*ones(n-1,1);
  56.     y = y*ones(n,1);
  57.   else
  58.     nx = x.nr;
  59.     ny = y.nr;
  60.     nz = z.nr;
  61.     if ((ny - nx - 1) || (ny - nz - 1)) {
  62.       error("Dimensions of vector arguments are incorrect.")
  63.     }
  64.   }
  65.  
  66.   T = diag(x, -1) + diag(y) + diag(z, 1);    // For non-sparse matrix.
  67.  
  68.   // n = max(size(y));                // For sparse matrix
  69.   // T = spdiags([ [x;0] y [0;z] ], -1:1, n, n);
  70.  
  71.   return T;
  72. };
  73.