home *** CD-ROM | disk | FTP | other *** search
/ Turbo Toolbox / Turbo_Toolbox.iso / 1988 / 02 / nonlin.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-11-23  |  5.5 KB  |  249 lines
  1.  
  2. (*---------------------------------------------------------------*)
  3. (*                            NONLIN.PAS                         *)
  4. (*          Loesung eines nichtlinearen Gleichungssystems        *)
  5. (*        (C) Prof.Dr.R.Rettig-Zinke & PASCAL INTERNATIONAL      *)
  6. (*---------------------------------------------------------------*)
  7. PROGRAM Nonlin;
  8.  
  9. CONST max = 10;                (*   max >= n + 1  !!!   *)
  10.  
  11. TYPE vektor = ARRAY[1..max] OF REAL;
  12.      matrix = ARRAY[1..max,1..max] OF REAL;
  13.  
  14. VAR   x, f, dx                : vektor;
  15.       a                       : matrix;
  16.       i, j, k,
  17.       m, m1, n,
  18.       iter, maxiter           : INTEGER;
  19.       schranke, sdiff, sneu   : REAL;
  20.       abl                     : BOOLEAN;
  21.  
  22. PROCEDURE Funktionen;
  23. (****************************************************************
  24.   hier: Gleichungssystem eingeben !!!
  25.  ****************************************************************)
  26. BEGIN
  27.   f[1] :=           x[1]+x[2]   -   x[1]*x[3]  +  0.5 *x[3] - 1.5;
  28.   f[2] :=           x[1]*x[2]   +   x[1]*x[1]  -  x[3]*x[3] +   6;
  29.   f[3] := Sqrt(x[1]*x[1]*x[1])  +   x[2]*x[2]  +  x[3]*x[3] -  14;
  30. END;
  31.  
  32.  
  33. PROCEDURE Ableitungen;
  34.  
  35. VAR xh,fh : vektor;
  36.     dj    : REAL;
  37.  
  38. BEGIN
  39. (***************************************************************
  40.   hier: evtl.  1. Ableitungen als Differentialquotienten
  41.                eingeben !!!
  42.  ***************************************************************)
  43.   IF abl THEN BEGIN
  44.     a[1,1] :=             1-x[3];
  45.     a[1,2] :=                  1;
  46.     a[1,3] :=          -x[1]+0.5;
  47.     a[2,1] :=        x[2]+2*x[1];
  48.     a[2,2] :=               x[1];
  49.     a[2,3] :=            -2*x[3];
  50.     a[3,1] :=     1.5*Sqrt(x[1]);
  51.     a[3,2] :=             2*x[2];
  52.     a[3,3] :=             2*x[3];
  53.   END
  54. (***************************************************************
  55.                1. Ableitungen als Differenzenquotienten
  56.                 (werden IN der Prozedur berechnet !!!)
  57.  ***************************************************************)
  58.   ELSE BEGIN
  59.     FOR j := 1 TO m DO
  60.       xh[j] := x[j];
  61.     Funktionen;
  62.     FOR i := 1 TO n DO
  63.       fh[i] := f[i];
  64.     FOR j := 1 TO m DO BEGIN
  65.       dj := ABS(x[j] / 1E4) + 1E-08;
  66.       x[j] := x[j] + dj;
  67.       Funktionen;
  68.       FOR i := 1 TO n DO
  69.         a[i,j] := f[i];
  70.       x[j] := xh[j];
  71.     END;
  72.     FOR j := 1 TO m DO BEGIN
  73.        dj := ABS(x[j] / 1E4) + 1E-08;
  74.        x[j] := x[j] - dj;
  75.        Funktionen;
  76.        FOR i := 1 TO n DO
  77.          a[i,j] := (a[i,j] - f[i]) / 2 / dj;
  78.        x[j] := xh[j];
  79.     END;
  80.     FOR j := 1 TO m DO
  81.       x[j] := xh[j];
  82.     FOR i := 1 TO n DO
  83.       f[i] := fh[i];
  84.   END;
  85. END;
  86.  
  87.  
  88. PROCEDURE Householder;
  89.  
  90. VAR   d         : vektor;
  91.       sigma, s,
  92.       sum, beta : REAL;
  93.       ok        : BOOLEAN;
  94.  
  95. BEGIN
  96.   ok := FALSE;
  97.   j := 1;
  98.   WHILE (j <= m) AND NOT (ok) DO BEGIN
  99.     sigma := 0;
  100.     FOR i := j TO n DO
  101.       sigma := sigma + Sqr( a[i,j] );
  102.     IF sigma <> 0 THEN BEGIN
  103.       IF a[j,j] < 0 THEN BEGIN
  104.         s := Sqrt(sigma);
  105.         d[j] := s;
  106.       END
  107.       ELSE BEGIN
  108.         s := - Sqrt(sigma);
  109.         d[j] := s;
  110.       END;
  111.       beta := 1 / (s * a[j,j] - sigma);
  112.       a[j,j] := a[j,j] - s;
  113.       FOR k := j+1 TO m1 DO BEGIN
  114.         sum := 0;
  115.         FOR i := j TO n DO
  116.           sum := sum + a[i,j] * a[i,k];
  117.         sum := beta * sum;
  118.         FOR i := j TO n DO
  119.           a[i,k] := a[i,k] + a[i,j] * sum;
  120.       END;
  121.     END
  122.     ELSE BEGIN
  123.       WriteLn;
  124.       WriteLn ('Singulaere Matrix !');
  125.       WriteLn;
  126.       ok := TRUE;
  127.       Halt;
  128.     END;
  129.     j := Succ(j);
  130.   END;
  131.   FOR j := m DOWNTO 1 DO BEGIN
  132.     dx[j] := a[j,m1];
  133.     FOR i := m DOWNTO j+1 DO
  134.       dx[j] := dx[j] - a[j,i] * dx[i];
  135.     dx[j] := dx[j] / d[j];
  136.   END;
  137. END;
  138.  
  139.  
  140. PROCEDURE ErweitereMatrix;
  141. BEGIN
  142.   FOR i := 1 TO n DO
  143.     a[i,m1] := f[i];
  144. END;
  145.  
  146.  
  147. PROCEDURE NeueWerte;
  148.  
  149. VAR sum : REAL;
  150.  
  151. BEGIN
  152.   FOR j := 1 TO m DO
  153.     x[j] := x[j] - dx[j];
  154. END;
  155.  
  156.  
  157. PROCEDURE Konvergenz;
  158.  
  159. VAR salt : REAL;
  160.  
  161. BEGIN
  162.   salt := sneu;
  163.   sneu := 0;
  164.   FOR i := 1 TO n DO
  165.     FOR j := 1 TO m DO
  166.       sneu := sneu + Sqr(a[i,j]);
  167.   sneu  := Sqrt(sneu);
  168.   sdiff := ABS(sneu - salt);
  169. END;
  170.  
  171.  
  172. PROCEDURE Bild;
  173. BEGIN
  174.   ClrScr;
  175.   Write ('iter    Schranke');
  176.   FOR j := 1 TO m DO
  177.     Write ('        x(',j,')');
  178.   WriteLn;
  179.   WriteLn;
  180. END;
  181.  
  182.  
  183. PROCEDURE Ausdruck;
  184. BEGIN
  185.   Write (iter:4, sdiff:12:6);
  186.   FOR j := 1 TO m DO
  187.     Write (x[j]:12:6);
  188.   WriteLn;
  189. END;
  190.  
  191.  
  192. PROCEDURE Ergebnis;
  193. BEGIN
  194.   WriteLn;
  195.   Write ('Bitte beliebige Taste druecken ... ');
  196.   REPEAT UNTIL KeyPressed;
  197.   ClrScr;
  198.   WriteLn ('Anzahl der Iterationen = ', iter);
  199.   WriteLn;
  200.   WriteLn ('  j        x(j)     dx(j)');
  201.   WriteLn;
  202.   FOR j := 1 TO m DO
  203.     WriteLn (j:3, x[j]:12:6, dx[j]:12:6);
  204.   WriteLn; WriteLn;
  205.   WriteLn ('  i        f(i)');
  206.   WriteLn;
  207.   FOR i := 1 TO n DO
  208.     WriteLn (i:3, f[i]:12:6);
  209.   WriteLn;
  210.   WriteLn ('Letzte Differenz: ',sdiff:12:6);
  211.   WriteLn;
  212. END;
  213.  
  214.  
  215. PROCEDURE Anfangswerte;
  216. BEGIN
  217.   n    :=  3;
  218.   m    :=  3;
  219.   x[1] := 15;
  220.   x[2] := 10;
  221.   x[3] :=  5;
  222.   maxiter  := 20;
  223.   schranke := 1E-6;
  224.   abl      := FALSE;
  225. END;
  226.  
  227. BEGIN (* Hauptprogramm *)
  228.   Anfangswerte;
  229.   m1    := m + 1;
  230.   iter  := 0;
  231.   sneu  := 0;
  232.   sdiff := schranke;
  233.   Bild;
  234.   Ausdruck;
  235.   WriteLn;
  236.   REPEAT
  237.     iter := Succ(iter);
  238.     Funktionen;
  239.     Ableitungen;
  240.     Konvergenz;
  241.     ErweitereMatrix;
  242.     Householder;
  243.     Ausdruck;
  244.     NeueWerte;
  245.   UNTIL (iter  =  maxiter)  OR
  246.         (sdiff < schranke);
  247.   Ergebnis;
  248. END.
  249.