home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / TURBOPAS / PAS_ENG.ZIP / GAUSID.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-07-19  |  3.2 KB  |  178 lines
  1. program gausid;  { -> 129 }
  2. { pascal program to perform simultaneous solution }
  3. { by Gauss-Seidel }
  4. { procedure SEID is included }
  5.  
  6. const maxr = 8;
  7.  maxc = 8;
  8.  
  9. type ary = array[1..maxr] of real;
  10.  arys = array[1..maxc] of real;
  11.  ary2s = array[1..maxr,1..maxc] of real;
  12.  
  13. var y : ary;
  14.  coef : arys;
  15.  a : ary2s;
  16.  n,m : integer;
  17.  first,
  18.  error : boolean;
  19.  
  20.  
  21. procedure get_data
  22.  (var a : ary2s;
  23.   var y : ary;
  24.   var n,m: integer);
  25. { get values for n and arrays a,y }
  26.  
  27. var i,j : integer;
  28.  
  29. begin
  30.   writeln;
  31.   repeat
  32.     write('How many equations? ');
  33.     readln(n);
  34.     if first then first:=false else clrscr
  35.   until n<maxr;
  36.   m:=n;
  37.   if n>1 then
  38.     begin
  39.       for i:=1 to n do
  40.  begin
  41.    writeln('Equation',i:3);
  42.    for j:=1 to n do
  43.      begin
  44.        write(j:3,':');
  45.        read(a[i,j])
  46.      end;
  47.  write(' C:');
  48.  read(y[i]);
  49.  readln  { clear the line }
  50.       end;
  51.     writeln;
  52.     for i:=1 to n do
  53.       begin
  54.  for j:=1 to m do
  55.    write(a[i,j]:7:4,' ');
  56.  writeln(':',y[i]:7:4)
  57.       end;
  58.       writeln
  59.     end  { if n>1 }
  60.   else if n<0 then n:=-n;
  61.   m:=n
  62. end;  { procedure get_data }
  63.  
  64. procedure write_data;
  65. { print out the answers }
  66.  
  67. var i : integer;
  68.  
  69. begin
  70.   for i:=1 to m do
  71.     write(coef[i]:9:5);
  72.   writeln
  73. end;  { write_data }
  74.  
  75. procedure seid
  76.  (a : ary2s;
  77.  y : ary;
  78.  var coef: arys;
  79.  ncol : integer;
  80.        var error: boolean);
  81. { matrix solution by Gauss Seidel }
  82.  
  83. const tol = 1.0E-4;
  84.  max = 100;
  85.  
  86. var done : boolean;
  87.        i,j,k,l,n: integer;
  88.  
  89.  nextc,hold,
  90.  sum,lambda,
  91.  ab,big : real;
  92.  
  93. begin
  94.   repeat
  95.     write('Relaxation factor? ');
  96.     readln(lambda)
  97.   until (lambda<2) and (lambda>0.0);
  98.   error:=false;
  99.   n:=ncol;
  100.   for i:=1 to n-1 do
  101.     begin
  102.       big:=abs(a[i,i]);
  103.       l:=i;
  104.       for j:=i+1 to n do
  105.  begin
  106.  { search for largest element }
  107.    ab:=abs(a[j,i]);
  108.    if ab>big then
  109.      begin
  110.        big:=ab;
  111.        l:=j
  112.      end
  113.       end;  { j-loop }
  114.     if big=0.0 then error:=true
  115.     else
  116.       begin
  117.  if l<>i then
  118.    begin
  119.    { interchange rows to put }
  120.    { largest element on diagonal }
  121.    for j:=1 to n do
  122.      begin
  123.        hold:=a[l,j];
  124.        a[l,j]:=a[i,j];
  125.        a[i,j]:=hold
  126.      end;
  127.    hold:=y[l];
  128.    y[l]:=y[i];
  129.    y[i]:=hold
  130.  end { if l<>i }
  131.       end { if big }
  132.     end; { i-loop }
  133.   if a[n,n]=0.0 then error:=true
  134.   else
  135.     begin
  136.       for i:=1 to n do
  137.  coef[i]:=0.0;  { initial guess }
  138.       i:=0;
  139.       repeat
  140.  i:=i+1;
  141.  done:=true;
  142.  for j:=1 to n do
  143.    begin
  144.      sum:=y[j];
  145.      for k:=1 to n do
  146.        if j<>k then
  147.   sum:=sum-a[j,k]*coef[k];
  148.      nextc:=sum/a[j,j];
  149.      if abs(nextc-coef[j])>tol then
  150.        begin
  151.   done:=false;
  152.   if nextc*coef[j]<0.0 then
  153.     nextc:=(coef[j]+nextc)*0.5
  154.        end;
  155.      coef[j]:=lambda*nextc+(1.0-lambda)*coef[j];
  156.      writeln(i:4,',coef(',j,')=',coef[j])
  157.  end  { j-loop }
  158.       until done or (i>max)
  159.   end; { if a[n,n]=0 }
  160.   if i>max then error:=true;
  161.   if error then writeln('ERROR: Matrix is singular')
  162. end;  { SEID }
  163.  
  164. begin  { MAIN program }
  165.   first:=true;
  166.   clrscr;
  167.   writeln;
  168.   writeln('Simultaneous solution by Gauss-Seidel');
  169.   repeat
  170.     get_data(a,y,n,m);
  171.     if n>1 then
  172.       begin
  173.  seid(a,y,coef,n,error);
  174.  if not error then write_data
  175.       end
  176.   until n<2
  177. end.
  178.