home *** CD-ROM | disk | FTP | other *** search
/ HAM Radio 1 / HamRadio.cdr / tech / eepub05 / ludec.pas < prev    next >
Pascal/Delphi Source File  |  1987-02-22  |  3KB  |  106 lines
  1. program ludec;     { LU Decomposition in Place }
  2. {                    Jeff Crawford     TRW     3-21-86
  3.  
  4. References:   "Computer Methods for Circuit Analysis and Design", Vlach
  5.               Singhal;  Van Nostrand Reinhold Company, 1983.
  6.  
  7.               "Applied Numerical Analysis", Curtis F. Gerald, Patrick
  8.               O. Wheatley, Addison-Wesley Publishing Co., 1984.
  9.  
  10. This program is unique compared to other LU Decompositions seen in the
  11. references above.  Due to the nature of the algorithm, it is possible to
  12. do the reduction "in-place" rather than identifying two additional matrices
  13. L & U, thus using more memory.
  14.  
  15. This program will readily allow computations on complex systems by
  16. structuring the input matrix as shown below:
  17.  
  18.           ! a11 a12  -b11 -b12 !  ! c11 c12 !  =  ! x11 x12 !
  19.           ! a21 a22  -b21 -b22 !  ! c21 c22 !  =  ! x21 x22 !
  20.           ! b11 b12   a11  a12 !  ! d11 d12 !  =  ! y11 y12 !
  21.           ! b21 b22   a21  a22 !  ! d21 d22 !  =  ! y21 y22 !
  22.  
  23. where in reality the equations are of the form
  24.  
  25. ! a11 +jb11  a12 + jb12 ! ! c11 + jd11  c12 + jd12 ! = ! x11 +jy11  x12 +jy12 !
  26. ! a21 +jb21  a22 + jb22 ! ! c21 + jd21  c22 + jd22 !   ! x21 +jy21  x22 +jy22 !
  27.                                                                               }
  28. type     mat1 = array[1..20,1..20] of real;
  29.  
  30. var      A         : mat1;
  31.          i,j,m,k   : integer;
  32.          n         : integer;
  33.          Sum       : real;
  34.          B,X       : array[1..10] of real;
  35.  
  36. begin
  37.  
  38. clrscr;
  39. textcolor(10);
  40. gotoxy(10,10);
  41. write('Enter Order of System  ');
  42. read(n);
  43. for i:= 1 to n do
  44.     begin
  45.     for j:= 1 to n do
  46.          begin
  47.          gotoxy(8,10);
  48.          clreol;
  49.          write('A[',i:2,' ',j:2,' ] = ');
  50.          read(A[i,j]);
  51.          end;
  52.     gotoxy(8,10);
  53.     clreol;
  54.     write('         B[',i,'] = ');
  55.     read(B[i]);
  56.     end;
  57. clrscr;
  58. for j:= 2 to n do  A[1,j]:= A[1,j]/A[1,1];
  59. for i:= 2 to n do                      { I is Main Outer Loop Counter }
  60.     begin
  61.     for k:= i to n do                  { K is Which Row of L Working on }
  62.          begin
  63.          Sum:= 0.0;
  64.          for j:= 1 to i-1 do           { J is # of Terms to Sum }
  65.               begin
  66.               Sum:= Sum + A[k,j]*A[j,i];
  67.               end;
  68.          A[k,i]:= A[k,i] - Sum;
  69.     end;
  70.  
  71.     for m:= i+1 to n do                { M = Which Column Row Element is in }
  72.          begin
  73.          Sum:= 0.0;
  74.          for j:= 1 to i-1 do  Sum:= Sum + A[i,j]*A[j,m];
  75.          A[i,m]:= (A[i,m] - Sum) / A[i,i];
  76.          end;
  77.     end;
  78.  
  79. B[1]:= B[1]/A[1,1];                    { Forward Substitution }
  80. for i:= 2 to n do
  81.     begin
  82.     Sum:= 0.0;
  83.     for k:= 1 to i-1 do      Sum:= Sum + A[i,k]*B[k];
  84.     B[i]:= ( B[i] - Sum ) / A[i,i];
  85.     end;
  86. X[n]:= B[n];
  87. j:= n-1;
  88. while j >= 0 do
  89.     begin
  90.     sum:= 0.0;
  91.     for k:= j+1 to n do      Sum:= Sum + A[j,k]*X[k];
  92.     X[j]:= B[j] - Sum;
  93.     j:= j - 1;
  94.     end;
  95. clrscr;
  96. gotoxy(1,4);
  97. writeln('                           Solution Vector   ');
  98. writeln('                       LU Decomposition Method ');
  99. writeln;
  100. for i:= 1 to n do writeln('                         X[ ',i,' ] = ',X[i]:10:6);
  101. end.
  102.  
  103.  
  104.  
  105.  
  106.