home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / wilk_r < prev   
Encoding:
Text File  |  1994-12-21  |  2.1 KB  |  74 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Various specific matrices devised/discussed by Wilkinson.
  4.  
  5. // Syntax:      wilk ( N )
  6.  
  7. // Description:
  8.  
  9. //      L = WILK (N) is the matrix or system of order N.
  10. //      N = 3: upper triangular system Ux=b illustrating inaccurate solution.
  11. //      N = 4: lower triangular system Lx=b, ill-conditioned.
  12. //      N = 5: HILB(6)(1:5,2:6)*1.8144.  Symmetric positive definite.
  13. //      N = 21: W21+, tridiagonal.   Eigenvalue problem.
  14.  
  15. //      References:
  16. //       J.H. Wilkinson, Error analysis of direct methods of matrix inversion,
  17. //          J. Assoc. Comput. Mach., 8 (1961),  pp. 281-330.
  18. //       J.H. Wilkinson, Rounding Errors in Algebraic Processes, Notes on Applied
  19. //          Science No. 32, Her Majesty's Stationery Office, London, 1963.
  20. //       J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University
  21. //          Press, 1965.
  22.  
  23. //    This file is a translation of wilk.m from version 2.0 of
  24. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  25. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  26.  
  27. // Dependencies
  28.    require hilb
  29.  
  30. //-------------------------------------------------------------------//
  31.  
  32. wilk = function (n)
  33. {
  34.   if (n == 3)
  35.   {
  36.     // Wilkinson (1961) p.323.
  37.     A = [ 1e-10,   .9,  -.4;
  38.            0,      .9,  -.4;
  39.            0,       0,  1e-10];
  40.     b = [  0,       0,    1]';
  41.     return << A=A; b=b >>;
  42.     
  43.   else if (n == 4) {
  44.  
  45.     // Wilkinson (1963) p.105.
  46.     A = [0.9143e-4,  0,          0,          0;
  47.           0.8762,    0.7156e-4,  0,          0;
  48.           0.7943,    0.8143,     0.9504e-4,  0;
  49.           0.8017,    0.6123,     0.7165,     0.7123e-4];
  50.     b = [0.6524,     0.3127,     0.4186,     0.7853]';
  51.     return << A=A; b=b >>;
  52.  
  53.   else if (n == 5) {
  54.  
  55.     // Wilkinson (1965), p.234.
  56.     A = hilb(6);
  57.     A = A[1:5; 2:6]*1.8144;
  58.     return A;
  59.  
  60.   else if (n == 21) {
  61.  
  62.     // Taken from gallery.m.  Wilkinson (1965), p.308.
  63.     E = diag(ones(n-1,1),1);
  64.     m = (n-1)/2;
  65.     A = diag(abs(-m:m)) + E + E';
  66.     return A;
  67.  
  68.   else
  69.  
  70.     error ("Sorry, that value of N is not available.");
  71.  
  72.   } } } }
  73. };
  74.