home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / LINEARAL.PAK / TRIDIAGO.M < prev   
Encoding:
Text File  |  1992-07-29  |  2.5 KB  |  96 lines
  1.  
  2.  
  3. (*:Version: Mathematica 2.0 *)
  4.  
  5. (*:Name: LinearAlgebra`Tridiagonal` *)
  6.  
  7. (*:Title: Tridiagonal Matrix Routines *)
  8.  
  9. (*:Author:
  10.     Jerry B. Keiper (Wolfram Research Inc.)
  11. *)
  12.  
  13. (*:Keywords:
  14.     tridiagonal matrix, diagonals
  15. *)
  16.  
  17. (*:Requirements: Length[a] >= Length[r] - 1,
  18.                  Length[b] >= Length[r],
  19.                  Length[c] >= Length[r] - 1. *)
  20.  
  21. (*:Warnings: Infinity introduced if pivot becomes zero. *)
  22.  
  23. (*:Sources:
  24.  
  25. *)
  26.  
  27. (*:Summary:
  28. TridiagonalSolve solves A . x == r for x, where A is an n x n
  29. tridiagonal matrix and r is a vector of length n.
  30. The 3 main diagonals of A are given by
  31.  
  32.      b[[1]] c[[1]]
  33.      a[[1]] b[[2]] c[[2]]
  34.             a[[2]]    .     .
  35.                       .     .     .
  36.  
  37.                                      c[[n-1]]
  38.                             a[[n-1]] b[[n]]  ,
  39.  
  40. Regular Gaussian elimination is employed, but no pivoting is done.
  41. *)
  42.  
  43. BeginPackage["LinearAlgebra`Tridiagonal`"]
  44.  
  45. TridiagonalSolve::usage =
  46. "TridiagonalSolve[ a, b, c, r ] solves A . x == r for x,
  47. where A is a tridiagonal matrix.  The three diagonals of A are given
  48. by: \n
  49.  
  50.      a11 = b[[1]], a12 = c[[1]],\n
  51.       a21 = a[[1]], a22 = b[[2]], a23 = c[[2]], ...\n
  52.  
  53. No pivoting is done, so the algorithm may fail even though a solution exists.
  54. This cannot happen for a symmetric positive definite matrix."
  55.  
  56. Begin["LinearAlgebra`Tridiagonal`Private`"]
  57.  
  58. TridiagonalSolve[a_List, b_List, c_List, r_List]:=
  59.     Module[{len=Length[r], solution=Array[0, Length[r]], aux, 
  60.             aux1=Array[0, Length[r]], a1 = Prepend[a, 0], iter}, 
  61.     aux = 1/b[[1]];
  62.     solution[[1]] = r[[1]] aux;
  63.     Do[aux1[[iter]] = c[[iter-1]] aux;
  64.         aux = 1/(b[[iter]]-a1[[iter]] aux1[[iter]]);
  65.         solution[[iter]] = (r[[iter]]-a1[[iter]] solution[[iter-1]]) aux,
  66.         {iter, 2, len}];             
  67.     Do[solution[[iter]] -= aux1[[iter+1]] solution[[iter+1]],
  68.         {iter, len-1, 1, -1}];
  69.     solution
  70.     ] /; Length[a] + 1 == Length[b] == Length[c] + 1 == Length[r]
  71.  
  72.  
  73. End[]  (* LinearAlgebra`Tridiagonal`Private` *)
  74.  
  75. Protect[ TridiagonalSolve ];
  76.  
  77. EndPackage[] (* LinearAlgebra`Tridiagonal` *)
  78.  
  79. (*:Limitations: No pivoting done. *)
  80.  
  81.  
  82. (*:Examples:
  83.  
  84. TridiagonalSolve[ Table[-1,{i,9}],Table[2,{i,10}],
  85.     Table[-1,{i,9}],{1,0,0,0,0,0,0,0,0,0}]
  86.  
  87. LinearSolve[Table[Switch[Abs[i-j], 0, 2, 1, -1, _, 0], {i, 10}, {j, 10}],
  88.     {1,0,0,0,0,0,0,0,0,0}]
  89.  
  90. TridiagonalSolve[{2,3,4,5}, {0,1,2,3,4}, {9,8,7,6}, {9,9,9,9,9}]
  91.  
  92. LinearSolve[ { { 0,9,0,0,0},{2,1,8,0,0},{0,3,2,7,0},{0,0,4,3,6},{0,0,0,5,4}},
  93.      {9,9,9,9,9}]
  94.  
  95. *)
  96.