home *** CD-ROM | disk | FTP | other *** search
/ Turbo Toolbox / Turbo_Toolbox.iso / turbo4 / hilb.pas < prev    next >
Pascal/Delphi Source File  |  1987-12-08  |  6KB  |  236 lines
  1. program Hilb;
  2. {
  3.  
  4.   Dieses Programm demonstriert die hohe Geschwindigkeit und
  5.   Genauigkeit der Real-Arithmetik von Turbo Pascal.
  6.  
  7.   --------------------------------------------------
  8.   Entnommen aus: Pascal Programs for Scientists and Engineers
  9.  
  10.   Alan R. Miller, Sybex
  11.   n x n Inverse Hilbert-Matrix
  12.   Lösung: 1 1 1 1 1
  13.   --------------------------------------------------
  14.  
  15.   Die Matrix wird über das Gauss-Jordansche Eliminationsverfahren
  16.   vereinfacht und simultan gelöst.
  17.  
  18.   1.  Setzen Sie den Compiler-Schalter O/C/Numeric processing auf "Software",
  19.       bevor Sie das Programm compilieren und starten.
  20.   2.  Wenn Ihr Computer mit einen numerischen Corprozessor ausgerüstet ist:
  21.       Compilieren Sie das Programm mit O/C/Numeric processing.. "Hardware"
  22.       ein zweites Mal und vergleichen Sie Geschwindigkeit und Genauigkeit.
  23. }
  24.  
  25. const
  26.   maxr = 10;
  27.   maxc = 10;
  28.  
  29. type
  30. {$IFOPT N+}                        { Wenn ein 80x87 installiert ist: }
  31.   real  = extended;                { Verwendung des Formats Extended }
  32. {$ENDIF}
  33.   ary   = array[1..maxr] of real;
  34.   arys  = array[1..maxc] of real;
  35.   ary2s = array[1..maxr, 1..maxc] of real;
  36.  
  37. var
  38.   y          : arys;
  39.   coef       : arys;
  40.   a, b       : ary2s;
  41.   n, m, i, j : integer;
  42.   error      : boolean;
  43.  
  44. procedure gaussj
  45.   (var b     : ary2s;  { Matrix mit den Quadraten der Koeffizienten }
  46.     y        : arys;  { Konstanten-Vektor (zur Inversion) }
  47.     var coef : arys;  { Lösungsvektor (Koeffizienten) }
  48.     ncol     : integer;  { Ordnung der Matrix (Spaltenzahl) }
  49.     var error: boolean); { TRUE, wenn die Matrix singulär ist }
  50.  
  51.  
  52. {  Gauss-Jordansche Inversion und Lösung, nach McCormick (8-FEB-81)
  53.    B(N,N) ist die Koeffizienten-Matrix und wird invertiert
  54.    Y(N)   ist der originale Konstante Vektor
  55.    W(N,M) besitzt konstante Vektoren und wird zur Lösung
  56.    Determ ist die Determinante der Matrix
  57.    Error  bekommt den Wert TRUE bei singulären Matrizen
  58.    nv     ist die Anzahl der konstanten Vektoren
  59. }
  60.  
  61. var
  62.   w    : array[1..maxc, 1..maxc] of real;
  63.   index: array[1..maxc, 1..3] of integer;
  64.   i, j, k, l, nv, irow, icol, n, l1   : integer;
  65.   determ, pivot, hold, sum, t, ab, big: real;
  66.  
  67. procedure swap(var a, b: real);
  68. var
  69.   hold: real;
  70. begin
  71.   hold := a; a := b; b := hold
  72. end; { swap }
  73.  
  74.  
  75. begin     { gaussj }
  76.   error := false;
  77.   nv := 1;   { Anzahl konstanter Vektoren }
  78.   n := ncol;
  79.   for i := 1 to n do
  80.     begin
  81.       w[i, 1] := y[i]; { Kopie des konstanten Vektors }
  82.       index[i, 3] := 0
  83.     end;
  84.   determ := 1.0;
  85.   for i := 1 to n do
  86.     begin
  87.       { Suche nach dem größten Element ("Drehpunkt") }
  88.       big := 0.0;
  89.       for j := 1 to n do
  90.           if index[j, 3] <> 1 then
  91.               for k := 1 to n do
  92.                 begin
  93.                   if index[k, 3] > 1 then
  94.                     begin
  95.                       writeln(' FEHLER: Matrix ist singulär');
  96.                       error := true;
  97.                       exit;
  98.                     end;
  99.                   if index[k, 3] < 1 then
  100.                     if abs(b[j, k]) > big then
  101.                       begin
  102.                         irow := j;
  103.                         icol := k;
  104.                         big := abs(b[j, k])
  105.                       end
  106.                 end; { for k }
  107.       index[icol, 3] := index[icol, 3] + 1;
  108.       index[i, 1] := irow;
  109.       index[i, 2] := icol;
  110.  
  111.   { Austausch der Reihen -> größtes Element auf die Diagonale }
  112.   if irow <> icol then
  113.     begin
  114.       determ := - determ;
  115.       for l := 1 to n do
  116.         swap(b[irow, l], b[icol, l]);
  117.       if nv > 0 then
  118.         for l := 1 to nv do
  119.           swap(w[irow, l], w[icol, l])
  120.     end; { if irow }
  121.  
  122.       { "Drehpunkt"-Reihe durch größtes Element dividieren }
  123.       pivot := b[icol, icol];
  124.       determ := determ * pivot;
  125.       b[icol, icol] := 1.0;
  126.       for l := 1 to n do
  127.         b[icol, l] := b[icol, l] / pivot;
  128.       if nv > 0 then
  129.         for l := 1 to nv do
  130.           w[icol, l] := w[icol, l] / pivot;
  131.       {  Reduktion der restlichen Reihen }
  132.       for l1 := 1 to n do
  133.         begin
  134.           if l1 <> icol then
  135.             begin
  136.               t := b[l1, icol];
  137.               b[l1, icol] := 0.0;
  138.               for l := 1 to n do
  139.                 b[l1, l] := b[l1, l] - b[icol, l] * t;
  140.               if nv > 0 then
  141.                 for l := 1 to nv do
  142.                   w[l1, l] := w[l1, l] - w[icol, l] * t;
  143.             end   { if l1 }
  144.         end
  145.     end { i loop };
  146.  
  147.   if error then exit;
  148.   { Austausch der Spalten }
  149.   for i := 1 to n do
  150.     begin
  151.       l := n - i + 1;
  152.       if index[l, 1] <> index[l, 2] then
  153.         begin
  154.           irow := index[l, 1];
  155.           icol := index[l, 2];
  156.           for k := 1 to n do
  157.             swap(b[k, irow], b[k, icol])
  158.         end { if index }
  159.     end  { i loop };
  160.   for k := 1 to n do
  161.     if index[k, 3] <> 1 then
  162.       begin
  163.         writeln(' FEHLER: Matrix ist singulär');
  164.         error := true;
  165.         exit;
  166.       end;
  167.   for i := 1 to n do
  168.     coef[i] := w[i, 1];
  169. end; { gaussj }
  170.  
  171.  
  172. procedure get_data(var a : ary2s;
  173.                    var y : arys;
  174.                    var n, m : integer);
  175.  
  176. { erzeugt die Originalmatrix }
  177.  
  178. var
  179.   i, j : integer;
  180.  
  181. begin
  182.   for i := 1 to n do
  183.     begin
  184.       a[n,i] := 1.0/(n + i - 1);
  185.       a[i,n] := a[n,i]
  186.     end;
  187.   a[n,n] := 1.0/(2*n -1);
  188.   for i := 1 to n do
  189.     begin
  190.       y[i] := 0.0;
  191.       for j := 1 to n do
  192.         y[i] := y[i] + a[i,j]
  193.     end;
  194.   writeln;
  195.   if n < 7 then
  196.     begin
  197.       for i:= 1 to n  do
  198.         begin
  199.           for j:= 1 to m do
  200.             write( a[i,j] :7:5, '  ');
  201.           writeln( ' : ', y[i] :7:5)
  202.         end;
  203.       writeln
  204.     end  { if n<7 }
  205. end;
  206.  
  207. procedure write_data;
  208.  
  209. { Ausgabe der Lösungen }
  210.  
  211. var
  212.   i : integer;
  213.  
  214. begin
  215.   for i := 1 to m do
  216.     write( coef[i] :13:9);
  217.   writeln;
  218. end;
  219.  
  220.  
  221. begin  { Hauptprogramm }
  222.   a[1,1] := 1.0;
  223.   n := 2;
  224.   m := n;
  225.   repeat
  226.     get_data (a, y, n, m);
  227.     for i := 1 to n do
  228.       for j := 1 to n do
  229.         b[i,j] := a[i,j] { Arbeits-Array initialisieren };
  230.     gaussj (b, y, coef, n, error);
  231.     if not error then write_data;
  232.     n := n+1;
  233.     m := n
  234.   until n > maxr;
  235. end.
  236.