home *** CD-ROM | disk | FTP | other *** search
/ Turbo Toolbox / Turbo_Toolbox.iso / spezial / 15 / numint / numint.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1989-05-19  |  9.9 KB  |  282 lines
  1. (* ------------------------------------------------------------------------- *)
  2. (*                                 NUMINT.PAS                                *)
  3. (*                 Unterprogramme zur numerischen Integration                *)
  4. (*                                                                           *)
  5. (*           Copyright (c) 1989   TOOLBOX  &  Karsten Gieselmann             *)
  6. (* ------------------------------------------------------------------------- *)
  7.  
  8. UNIT NumInt;
  9.  
  10. INTERFACE
  11.  
  12.   CONST
  13.     MaxDim  = 3;                (* Maximale Dimension eines Mehrfachintegrals *)
  14.  
  15.   TYPE
  16. {$IFDEF CPU87}
  17.     FLOAT  = DOUBLE;                                        (* Fließkommatyp *)
  18. {$ELSE}
  19.     FLOAT  = REAL;
  20. {$ENDIF}
  21.     VECTOR = ARRAY[1..MaxDim] OF FLOAT;
  22.     POINTS = ARRAY[1..MaxDim] OF WORD;                       (* Stützstellen *)
  23.     PROC   = POINTER;{ TO FUNCTION(FLOAT ) : FLOAT }    (* Prozedurparameter *)
  24.     VPROC  = POINTER;{ TO FUNCTION(VECTOR) : FLOAT }
  25.  
  26.  
  27.   (* ------------- Grundlegende Integrationsverfahren (Regeln) ------------- *)
  28.  
  29.   FUNCTION Trapezoid(Integrand : PROC; a,b : FLOAT; N : WORD) : FLOAT;
  30.     (* Integration mit der Trapez-Regel *)
  31.  
  32.   FUNCTION Simpson(Integrand : PROC; a,b : FLOAT; N : WORD) : FLOAT;
  33.     (* Integration mit der Simpson-Regel *)
  34.  
  35.   FUNCTION Gauss(Integrand : PROC; a,b : FLOAT; N : WORD) : FLOAT;
  36.     (* Integration durch Gauss-Quadratur *)
  37.  
  38.  
  39.   (* --- Aufbauende Methoden zur automatischen Berechnung eines Integrals -- *)
  40.  
  41.   FUNCTION Successive(    Rule         : PROC;          (* Integrationsregel *)
  42.                           f            : PROC;        (* Integrandenfunktion *)
  43.                           a,b          : FLOAT;       (* Integrationsgrenzen *)
  44.                           Accuracy     : FLOAT;      (* relative Genauigkeit *)
  45.                           MaxPoints    : WORD;  (* Maximale Stützstellenzahl *)
  46.                       VAR Result       : FLOAT;        (* Wert des Integrals *)
  47.                       VAR Error        : FLOAT)          (* absoluter Fehler *)
  48.                                        : BOOLEAN;   (* Genauigkeit erreicht? *)
  49.  
  50.   FUNCTION Adaptive  (    Rule         : PROC;          (* Integrationsregel *)
  51.                           f            : PROC;        (* Integrandenfunktion *)
  52.                           a,b          : FLOAT;       (* Integrationsgrenzen *)
  53.                           Accuracy     : FLOAT;      (* relative Genauigkeit *)
  54.                           Points       : WORD;           (* Stützstellenzahl *)
  55.                           MaxIntervals : WORD;     (* maximale Intervallzahl *)
  56.                       VAR Result       : FLOAT;        (* Wert des Integrals *)
  57.                       VAR Error        : FLOAT)          (* absoluter Fehler *)
  58.                                        : BOOLEAN;   (* Genauigkeit erreicht? *)
  59.  
  60.  
  61.   (* ---------------- Berechnung von Mehrfach-Integralen ------------------- *)
  62.  
  63.   FUNCTION Romberg(       Integrand    : PROC;        (* Integrandenfunktion *)
  64.                           Dim          : WORD;  (* Dimension des Integranden *)
  65.                           a,b          : VECTOR;      (* Integrationsgrenzen *)
  66.                           Accuracy     : FLOAT;      (* relative Genauigkeit *)
  67.                           N            : POINTS;         (* Gewichtsfaktoren *)
  68.                           MaxStep      : WORD;       (* Maximale Schrittzahl *)
  69.                       VAR Result       : FLOAT;        (* Wert des Integrals *)
  70.                       VAR Error        : FLOAT)          (* absoluter Fehler *)
  71.                                        : BOOLEAN;   (* Genauigkeit erreicht? *)
  72.  
  73.  
  74.  
  75. IMPLEMENTATION
  76.  
  77.  
  78.   FUNCTION Trapezoid;
  79.     (* Integration mit der Trapez-Regel *)
  80.   VAR
  81.     h,T : FLOAT;
  82.     i   : WORD;
  83.  
  84.   FUNCTION f(x : FLOAT) : FLOAT;
  85.     INLINE($FF/$9E/Integrand);                     (* CALL FAR Integrand[BP] *)
  86.  
  87.   BEGIN
  88.     h := (b-a)/N;
  89.     T := (f(a)+f(b)) * 0.5;
  90.     FOR i:=1 TO N-1 DO
  91.       T := T + f(a+i*h);
  92.     Trapezoid := T*h
  93.   END;
  94.  
  95.  
  96.   FUNCTION Simpson;
  97.     (* Integration mit der Simpson-Regel *)
  98.   VAR
  99.     h,S : FLOAT;
  100.     i   : WORD;
  101.  
  102.   FUNCTION f(x : FLOAT) : FLOAT;
  103.     INLINE($FF/$9E/Integrand);                     (* CALL FAR Integrand[BP] *)
  104.  
  105.   BEGIN
  106.     h := (b-a)/(2*N);
  107.     S := f(a) + f(b);
  108.     FOR i:=0 TO N-1 DO
  109.       S := S + 4.0*f(a+(2*i+1)*h);
  110.     FOR i:=1 TO N-1 DO
  111.       S := S + 2.0*f(a+(2*i)*h);
  112.     Simpson := S * h/3.0;
  113.   END;
  114.  
  115.  
  116.   FUNCTION Gauss;
  117.     (* Integration durch Gauss-Quadratur *)
  118.   CONST
  119.     z : ARRAY[1..16] OF FLOAT = (                            (* Stützstellen *)
  120.         0.095012509837637440185,   -0.095012509837637440185,
  121.         0.281603550779258913230,   -0.281603550779258913230,
  122.         0.458016777657227386342,   -0.458016777657227386342,
  123.         0.617876244402643748447,   -0.617876244402643748447,
  124.         0.755404408355003033895,   -0.755404408355003033895,
  125.         0.865631202387831743880,   -0.865631202387831743880,
  126.         0.944575023073232576078,   -0.944575023073232576078,
  127.         0.989400934991649932596,   -0.989400934991649932596);
  128.  
  129.     w : ARRAY[1..16] OF FLOAT = (                                (* Gewichte *)
  130.          0.189450610455068496285,   0.189450610455068496285,
  131.          0.182603415044923588867,   0.182603415044923588867,
  132.          0.169156519395002538189,   0.169156519395002538189,
  133.          0.149595988816576732081,   0.149595988816576732081,
  134.          0.124628971255533872052,   0.124628971255533872052,
  135.          0.095158511682492784810,   0.095158511682492784810,
  136.          0.062253523938647892863,   0.062253523938647892863,
  137.          0.027152459411754094852,   0.027152459411754094852);
  138.   VAR
  139.     i,j   : WORD;
  140.     x,g,h : FLOAT;
  141.  
  142.   FUNCTION f(x : FLOAT) : FLOAT;
  143.     INLINE($FF/$9E/Integrand);                     (* CALL FAR Integrand[BP] *)
  144.  
  145.   BEGIN
  146.     h := (b-a)/(2*N);
  147.     g := 0;
  148.     FOR i:=1 TO N DO
  149.       FOR j:=1 TO 16 DO
  150.         BEGIN
  151.           x := a + h*(2*i-1+z[j]);
  152.           g := g + f(x)*w[j];
  153.         END;
  154.     Gauss := g * h
  155.   END;
  156.  
  157.  
  158.   FUNCTION Successive;
  159.     (* Integration durch sukzessive Stützstellenverdopplung *)
  160.   VAR
  161.     LastResult : FLOAT;
  162.     N          : WORD;
  163.     Success    : BOOLEAN;
  164.  
  165.   FUNCTION Integrate(f : PROC; a,b : FLOAT; N : WORD) : FLOAT;
  166.     INLINE($FF/$9E/Rule);                               (* CALL FAR Rule[BP] *)
  167.  
  168.   BEGIN
  169.     N := 2;
  170.     Success := FALSE;
  171.     Result := Integrate(f, a, b, N);             (* erste Näherung bestimmen *)
  172.     REPEAT
  173.       LastResult := Result;
  174.       N := N * 2;                                (* Stützstellen verdoppeln! *)
  175.       Result := Integrate(f, a, b, N);
  176.       Error := ABS(Result-LastResult);
  177.       Success := (Error < ABS(Accuracy*Result)); (* bis Genauigkeit erreicht *)
  178.     UNTIL (N > MaxPoints) OR Success;
  179.     Successive := Success;
  180.   END;
  181.  
  182.  
  183.   FUNCTION Adaptive;
  184.     (* Adaptive Integration durch rekursive Intervallhalbierung *)
  185.   CONST
  186.     Rule_ : PROC = NIL;
  187.   VAR
  188.     r,e,EstimateError : FLOAT;
  189.     Intervals         : WORD;
  190.  
  191.   FUNCTION Integrate(f : PROC; a,b : FLOAT; N : WORD) : FLOAT;
  192.     INLINE($FF/$1E/Rule_);                             (* CALL FAR Rule_[DS] *)
  193.  
  194.     FUNCTION SubInterval(a, b : FLOAT) : BOOLEAN;
  195.       (* rekusive Integration, liefert TRUE bei erfolgreicher Bearbeitung *)
  196.     BEGIN
  197.       SubInterval := FALSE;
  198.       IF Intervals <= MaxIntervals THEN BEGIN
  199.         r := Integrate(f, a, b, Points*2);
  200.         e := ABS(r - Integrate(f, a, b, Points));
  201.         IF e <= EstimateError THEN BEGIN            (* Genauigkeit erreicht? *)
  202.           Result := Result + r;                (* Teilergebnisse aufaddieren *)
  203.           Error  := Error  + e;
  204.         END ELSE BEGIN
  205.           INC(Intervals);
  206.           IF NOT SubInterval(a, (a+b)/2) THEN EXIT;
  207.           IF NOT SubInterval((a+b)/2, b) THEN EXIT;
  208.         END;
  209.         SubInterval := TRUE;
  210.       END;
  211.     END;
  212.  
  213.   BEGIN
  214.     Rule_ := Rule;     (* Zugriff über DS, da BP sich verändert (Rekursion!) *)
  215.     Intervals := 0;
  216.     Result := 0;
  217.     Error := 0;
  218.     EstimateError := ABS(Accuracy * Integrate(f, a, b, Points));
  219.     Adaptive := SubInterval(a, b);
  220.   END;
  221.  
  222.  
  223.   FUNCTION Romberg;
  224.     (* Berechnung eines Mehrfachintegrals nach dem Romberg-Verfahren *)
  225.   CONST
  226.     Integrand_   : PROC = NIL;
  227.   VAR
  228.     p,x          : VECTOR;
  229.     T            : ARRAY[1..100] OF FLOAT;
  230.     i,j,k        : WORD;
  231.     LastResult,w : FLOAT;
  232.  
  233.     FUNCTION f(x : VECTOR) : FLOAT;
  234.       INLINE($FF/$1E/Integrand_);                 (* CALL FAR Integrand_[DS] *)
  235.  
  236.     FUNCTION Sum(k : WORD) : FLOAT;
  237.       (* Rekursives Aufsummieren *)
  238.     VAR
  239.       i : WORD;
  240.       s : FLOAT;
  241.     BEGIN
  242.       s := 0;
  243.       FOR i:=1 TO N[k] DO BEGIN
  244.         x[k] := a[k] + p[k]*(i-0.5);
  245.         IF k > 1 THEN                          (* schon auf unterster Ebene? *)
  246.           s := s + Sum(k-1)
  247.         ELSE
  248.           s := s + f(x);
  249.         END;
  250.       Sum := s;
  251.     END;
  252.  
  253.   BEGIN
  254.     Romberg := FALSE;
  255.     IF NOT(Dim IN [1..MaxDim]) THEN EXIT;            (* unerlaubte Dimension *)
  256.     Integrand_ := Integrand;
  257.     LastResult := 0.0;
  258.     k := 1;
  259.     REPEAT
  260.       FOR i:=1 TO Dim DO BEGIN
  261.         p[i] := (b[i]-a[i]) / N[i];               (* Bestimmung der Gewichte *)
  262.       END;
  263.       T[k] := Sum(Dim);                   (* Aufsummieren der Funktionswerte *)
  264.       FOR i:=1 TO Dim DO
  265.         T[k] := T[k] * p[i];             (* Multiplikation mit den Gewichten *)
  266.       FOR j:=2 TO k DO BEGIN
  267.         i := Succ(k) - j;
  268.         w := 2 shl (j-2);
  269.         T[i] := T[i+1] + (T[i+1]-T[i])/(Sqr(w)-1.0);
  270.       END;
  271.       Result := T[1];
  272.       Error := ABS(Result - LastResult);            (* erreichte Genauigkeit *)
  273.       INC(k);
  274.       FOR i:=1 TO Dim DO   (* Zahl der Stützstellen pro Dimension verdoppeln *)
  275.         N[i] := 2*N[i];
  276.       LastResult := Result;
  277.     UNTIL (k > MaxStep) OR ((k > 2) AND (Error < ABS(Result*Accuracy)));
  278.     Romberg := Error < ABS(Result*Accuracy);
  279.   END;
  280.  
  281. END.
  282.