home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 1.ddi / CHAP5.ARC / INTEGRAT.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  41.5 KB  |  985 lines
  1. unit Integrat;
  2.  
  3. {----------------------------------------------------------------------------}
  4. {-                                                                          -}
  5. {-     Turbo Pascal Numerical Methods Toolbox                               -}
  6. {-     Copyright (c) 1986, 87 by Borland International, Inc.                -}
  7. {-                                                                          -}
  8. {-  This unit provides procedures for performing numerical integration.     -}
  9. {-                                                                          -}
  10. {----------------------------------------------------------------------------}
  11.  
  12. {$I Float.inc} { Determines the setting of the $N compiler directive }
  13.  
  14. interface
  15.  
  16. {$IFOPT N+}
  17. type
  18.   Float = Double; { 8 byte real, requires 8087 math chip }
  19. {$ELSE}
  20. type
  21.   Float = real;   { 6 byte real, no math chip required }
  22. {$ENDIF}
  23.  
  24. const
  25.   TNArraySize = 50;              { Size of the vectors }
  26.  
  27. type
  28.   TNvector = array[0..TNArraySize] of Float;
  29.  
  30. procedure Simpson(LowerLimit   : Float;
  31.                   UpperLimit   : Float;
  32.                   NumIntervals : integer;
  33.               var Integral     : Float;
  34.               var Error        : byte;
  35.                   FuncPtr      : Pointer);
  36.  
  37. {----------------------------------------------------------------------------}
  38. {-                                                                          -}
  39. {-  Purpose: given a function, TNTargetF(X), compute the integral of        -}
  40. {-           TNTargetF(X) from LowerLimit to UpperLimit using Simpson       -}
  41. {-           Composite Algorithm.                                           -}
  42. {-                                                                          -}
  43. {-  User-defined Functions: TNTargetF(X : real) : real;                     -}
  44. {-                                                                          -}
  45. {-    Global Variables:  LowerLimit   : real;    Lower limit of integration -}
  46. {-                       UpperLimit   : real;    Upper limit of integration -}
  47. {-                       NumIntervals : integer; Number of subintervals     -}
  48. {-                                               over which to use          -}
  49. {-                                               Simpson's rule.            -}
  50. {-                       Integral     : real;    Value of the integral of   -}
  51. {-                                               TNTargetF over the given   -}
  52. {-                                               interval                   -}
  53. {-                       Error        : byte;    Flags if something goes    -}
  54. {-                                               wrong                      -}
  55. {-                                                                          -}
  56. {-              Errors:  0: No errors                                       -}
  57. {-                       1: NumIntervals < 0                                -}
  58. {-                                                                          -}
  59. {----------------------------------------------------------------------------}
  60.  
  61. procedure Trapezoid(LowerLimit   : Float;
  62.                     UpperLimit   : Float;
  63.                     NumIntervals : integer;
  64.                 var Integral     : Float;
  65.                 var Error        : byte;
  66.                     FuncPtr      : Pointer);
  67.  
  68. {----------------------------------------------------------------------------}
  69. {-                                                                          -}
  70. {-  Purpose: Given a function, TNTargetF(X), compute the integral of        -}
  71. {-  TNTargetF(X) from LowerLimit to UpperLimit using the Trapezoid Rule     -}
  72. {-                                                                          -}
  73. {-  User-defined Functions: TNTargetF(X : real) : real;                     -}
  74. {-                                                                          -}
  75. {-    Global Variables:  LowerLimit   : real;    Lower limit of integration -}
  76. {-                       UpperLimit   : real;    Upper limit of integration -}
  77. {-                       NumIntervals : integer; Number of subintervals     -}
  78. {-                                               over which to use the      -}
  79. {-                                               trapezoid rule.            -}
  80. {-                       Integral     : real;    Value of the integral of   -}
  81. {-                                               TNTargetF over the given   -}
  82. {-                                               interval                   -}
  83. {-                       Error        : byte;    Flags if something goes    -}
  84. {-                                               wrong                      -}
  85. {-                                                                          -}
  86. {-              Errors:  0: No errors                                       -}
  87. {-                       1: NumIntervals < 0                                -}
  88. {-                                                                          -}
  89. {----------------------------------------------------------------------------}
  90.  
  91. procedure Adaptive_Simpson(LowerLimit   : Float;
  92.                            UpperLimit   : Float;
  93.                            Tolerance    : Float;
  94.                            MaxIntervals : integer;
  95.                        var Integral     : Float;
  96.                        var NumIntervals : integer;
  97.                        var Error        : byte;
  98.                            FuncPtr      : Pointer);
  99.  
  100. {----------------------------------------------------------------------------}
  101. {-                                                                          -}
  102. {-    Input: LowerLimit, UpperLimit, Tolerance, MaxIntervals                -}
  103. {-   Output: Integral, NumIntervals, Error                                  -}
  104. {-                                                                          -}
  105. {-  Purpose: Given a function, TNTargetF(X), compute the integral of        -}
  106. {-           TNTargetF from LowerLimit to UpperLimit using Adaptive         -}
  107. {-           Quadrature and Simpson's rule.                                 -}
  108. {-                                                                          -}
  109. {-  User-defined Functions: TNTargetF(X : real) : real;                     -}
  110. {-                                                                          -}
  111. {-  Global Variables:  LowerLimit   : real;    Lower limit of integration   -}
  112. {-                     UpperLimit   : real;    Upper limit of integration   -}
  113. {-                     Tolerance    : real;    Tolerance in answer          -}
  114. {-                     MaxIntervals : integer; Max. number of subintervals  -}
  115. {-                                             used to approximate integral -}
  116. {-                     Integral     : real;    Value of the integral of     -}
  117. {-                                             TNTargetF over the given     -}
  118. {-                                             interval                     -}
  119. {-                     NumIntervals : integer; Number of subintervals       -}
  120. {-                     Error        : byte;    Flags if something goes      -}
  121. {-                                             wrong                        -}
  122. {-                                                                          -}
  123. {-            Errors:  0: No errors                                         -}
  124. {-                     1: Tolerance <= 0                                    -}
  125. {-                     2: MaxIntervals <= 0                                 -}
  126. {-                     3: NumIntervals >= MaxIntervals                      -}
  127. {-                                                                          -}
  128. {-              Note:  Adaptive quadrature is a recursive routine.          -}
  129. {-                     In order to avoid recursive procedure calls,         -}
  130. {-                     (which slows down the execution) a stack is          -}
  131. {-                     created to simulate recursion.  The stack is         -}
  132. {-                     created using pointer implementation.                -}
  133. {-                                                                          -}
  134. {----------------------------------------------------------------------------}
  135.  
  136. procedure Adaptive_Gauss_Quadrature(LowerLimit   : Float;
  137.                                     UpperLimit   : Float;
  138.                                     Tolerance    : Float;
  139.                                     MaxIntervals : integer;
  140.                                 var Integral     : Float;
  141.                                 var NumIntervals : integer;
  142.                                 var Error        : byte;
  143.                                     FuncPtr      : Pointer);
  144.  
  145. {----------------------------------------------------------------------------}
  146. {-                                                                          -}
  147. {-  Purpose: Given a function, TNTargetF(X), compute the integral of        -}
  148. {-           TNTargetF from LowerLimit to UpperLimit using Adaptive         -}
  149. {-           Quadrature and Gaussian Quadrature with a 16th order           -}
  150. {-           Legendre polynomial.                                           -}
  151. {-                                                                          -}
  152. {-  User-defined Functions: TNTargetF(X : real) : real;                     -}
  153. {-                                                                          -}
  154. {-  Global Variables:  LowerLimit   : real;    Lower limit of integration   -}
  155. {-                     UpperLimit   : real;    Upper limit of integration   -}
  156. {-                     Tolerance    : real;    Tolerance in answer          -}
  157. {-                     MaxIntervals : integer; Max no. of subintervals over -}
  158. {-                                             Which to approximate integral-}
  159. {-                     Integral     : real;    Value of the integral of     -}
  160. {-                                             TNTargetF over the given     -}
  161. {-                                             interval                     -}
  162. {-                     NumIntervals : integer; Number of subintervals       -}
  163. {-                                             used to approximate integral -}
  164. {-                     Error        : byte;    Flags if something goes      -}
  165. {-                                             wrong                        -}
  166. {-                                                                          -}
  167. {-              Errors:  0: No errors                                       -}
  168. {-                       1: Tolerance <= 0                                  -}
  169. {-                       2: MaxIntervals <= 0                               -}
  170. {-                       3: NumIntervals >= MaxIntervals                    -}
  171. {-                                                                          -}
  172. {-                Note:  Adaptive quadrature is a recursive routine.        -}
  173. {-                       In order to avoid recursive procedure calls,       -}
  174. {-                       (which slow down the execution) a stack is         -}
  175. {-                       created to simulate recursion.  The stack is       -}
  176. {-                       created using pointer implementation.              -}
  177. {-                                                                          -}
  178. {----------------------------------------------------------------------------}
  179.  
  180. procedure Romberg(LowerLimit : Float;
  181.                   UpperLimit : Float;
  182.                   Tolerance  : Float;
  183.                   MaxIter    : integer;
  184.               var Integral   : Float;
  185.               var Iter       : integer;
  186.               var Error      : byte;
  187.                   FuncPtr    : Pointer);
  188.  
  189. {----------------------------------------------------------------------------}
  190. {-                                                                          -}
  191. {-    Input: LowerLimit, UpperLimit, Tolerance, MaxIter                     -}
  192. {-   Output: Integral, Iter, Error                                          -}
  193. {-                                                                          -}
  194. {-  Purpose: Given a function, TNTargetF(X), this procedure approximates    -}
  195. {-           the integral of TNTargetF from LowerLimit to UpperLimit        -}
  196. {-           using the Romberg method.                                      -}
  197. {-                                                                          -}
  198. {-  User-defined Functions: TNTargetF(X : real) : real;                     -}
  199. {-                                                                          -}
  200. {-    Global Variables:  LowerLimit : real;    Lower limit of integration   -}
  201. {-                       UpperLimit : real;    Upper limit of integration   -}
  202. {-                       Tolerance  : real;    Tolerance in answer          -}
  203. {-                       MaxIter    : integer; Maximum number of iterations -}
  204. {-                       Integral   : real;    Value of the integral of     -}
  205. {-                                             TNTargetF over the given     -}
  206. {-                                             interval                     -}
  207. {-                       Iter       : integer; Number of iterations         -}
  208. {-                       Error      : byte;    Flags if something goes      -}
  209. {-                                             wrong                        -}
  210. {-                                                                          -}
  211. {-              Errors:  0: No errors                                       -}
  212. {-                       1: Tolerance <= 0                                  -}
  213. {-                       2: MaxIter <= 0                                    -}
  214. {-                       3: Iter > MaxIter                                  -}
  215. {-                                                                          -}
  216. {-        Version Date:  26 January 1987                                    -}
  217. {-                                                                          -}
  218. {----------------------------------------------------------------------------}
  219.  
  220. implementation
  221.  
  222. {$F+}
  223. {$L Integrat.OBJ}
  224. function UserFunction(X : Float; ProcAddr : Pointer) : Float; external;
  225. {$F-}
  226.  
  227. procedure Simpson{(LowerLimit  : Float;
  228.                   UpperLimit   : Float;
  229.                   NumIntervals : integer;
  230.               var Integral     : Float;
  231.               var Error        : byte;
  232.                   FuncPtr      : Pointer)};
  233.  
  234. var
  235.   Spacing : Float;           { Size of each subinterval }
  236.   Point : Float;             { Midway point of each interval }
  237.   OddSum, EvenSum : Float;   { Sums of values over odd numbered }
  238.                             { and even numbered subintervals }
  239.   LimitsValue : Float;       { Sum of values at the endpoints }
  240.   Interval : integer;       { Counter }
  241.  
  242. begin { procedure Simpson }
  243.   if NumIntervals <= 0 then
  244.     Error := 1
  245.   else
  246.     begin
  247.       Spacing := (UpperLimit - LowerLimit) / (2 * NumIntervals);
  248.       Point := LowerLimit;
  249.       OddSum := 0;
  250.       EvenSum := 0;
  251.       for Interval := 1 to 2*NumIntervals - 1 do
  252.       begin
  253.         Point := Point + Spacing;
  254.         if Odd(Interval) then
  255.           OddSum := OddSum + UserFunction(Point, FuncPtr)
  256.         else
  257.           EvenSum := EvenSum + UserFunction(Point, FuncPtr);
  258.       end;
  259.       LimitsValue := UserFunction(UpperLimit, FuncPtr) +
  260.                      UserFunction(LowerLimit, FuncPtr);
  261.       Integral := Spacing * (LimitsValue + 2 * EvenSum + 4 * OddSum) / 3;
  262.     end;
  263. end; { procedure Simpson }
  264.  
  265. procedure Trapezoid{(LowerLimit  : Float;
  266.                     UpperLimit   : Float;
  267.                     NumIntervals : integer;
  268.                 var Integral     : Float;
  269.                 var Error        : byte;
  270.                     FuncPtr      : Pointer)};
  271.  
  272. var
  273.   Spacing : Float;           { Size of each subinterval }
  274.   Point : Float;             { Midway point of each interval }
  275.   Sum : Float;               { Sums of values over intervals }
  276.   LimitsValue : Float;       { Sum of values at the endpoints }
  277.   Interval : integer;       { Counter }
  278.  
  279. begin  { procedure Trapezoid }
  280.   if NumIntervals <= 0 then
  281.     Error := 1
  282.   else
  283.     begin
  284.       Spacing := (UpperLimit - LowerLimit) / NumIntervals;
  285.       Point := LowerLimit;
  286.       Sum := 0;
  287.       for Interval := 1 to NumIntervals - 1 do
  288.       begin
  289.         Point := Point + Spacing;
  290.         Sum := Sum + UserFunction(Point, FuncPtr);
  291.       end;
  292.       LimitsValue := UserFunction(UpperLimit, FuncPtr) +
  293.                      UserFunction(LowerLimit, FuncPtr);
  294.       Integral := Spacing * (LimitsValue + 2 * Sum) / 2;
  295.     end;
  296. end; { procedure Trapezoid }
  297.  
  298. procedure Adaptive_Simpson{(LowerLimit  : Float;
  299.                            UpperLimit   : Float;
  300.                            Tolerance    : Float;
  301.                            MaxIntervals : integer;
  302.                        var Integral     : Float;
  303.                        var NumIntervals : integer;
  304.                        var Error        : byte;
  305.                            FuncPtr      : Pointer)};
  306. type
  307.   IntegralData = record
  308.                    LB, VLB, UB, VUB, Est : Float;
  309.                  end;
  310.  
  311.   Ptr = ^StackItem;
  312.  
  313.   StackItem = record
  314.                 Info : IntegralData;
  315.                 Next : Ptr;
  316.               end;
  317.  
  318. var
  319.   Spacing, Estimate, NewEstimate1, NewEstimate2, NewEstimate,
  320.   Middle, ValueMiddle : Float;
  321.   Data : IntegralData;
  322.   Stack : Ptr;
  323.   Finished : boolean;
  324.  
  325. procedure InitializeStack(var Stack : Ptr);
  326. begin
  327.   Stack := nil;
  328. end; { procedure InitializeStack}
  329.  
  330. procedure Push(Data  : IntegralData;
  331.            var Stack : Ptr);
  332.  
  333. {------------------------------}
  334. {- Input: Data, Stack         -}
  335. {- Output: Stack              -}
  336. {-                            -}
  337. {- Push Data onto the Stack   -}
  338. {------------------------------}
  339.  
  340. var
  341.   NewNode : Ptr;
  342.  
  343. begin
  344.   New(NewNode);
  345.   NewNode^.Info := Data;
  346.   NewNode^.Next := Stack;
  347.   Stack := NewNode;
  348. end; { procedure Push }
  349.  
  350. procedure Pop(var Data  : IntegralData;
  351.               var Stack : Ptr);
  352.  
  353. {------------------------------}
  354. {- Input: Stack               -}
  355. {- Output: Data, Stack        -}
  356. {-                            -}
  357. {- Pop Data from the Stack    -}
  358. {------------------------------}
  359.  
  360. var
  361.   OldNode : Ptr;
  362.  
  363. begin
  364.   OldNode := Stack;
  365.   Stack := OldNode^.Next;
  366.   Data := OldNode^.Info;
  367.   Dispose(OldNode);
  368. end; { procedure Pop }
  369.  
  370. function Simpson(LowerLimit  : Float;
  371.                  ValueLowLmt : Float;
  372.                  UpperLimit  : Float;
  373.                  ValueUpLmt  : Float) : Float;
  374.  
  375. {----------------------------------------------------------}
  376. {- Input: LowerLimit, ValueLowLmt, UpperLimit, ValueUpLmt -}
  377. {-                                                        -}
  378. {- This function returns the integral of TNTargetF(X)     -}
  379. {- from LowerLimit to UpperLimit using Simpson's rule.    -}
  380. {----------------------------------------------------------}
  381.  
  382. var
  383.   Spacing : Float;
  384.  
  385. begin
  386.   Spacing := (UpperLimit - LowerLimit) / 2;
  387.   Simpson := Spacing * (ValueLowLmt + 4 * UserFunction(LowerLimit + Spacing, FuncPtr)
  388.              + ValueUpLmt) / 3;
  389. end; { function Simpson }
  390.  
  391. procedure TestAndInitialize(LowerLimit   : Float;
  392.                             UpperLimit   : Float;
  393.                             Tolerance    : Float;
  394.                             MaxIntervals : integer;
  395.                         var Data         : IntegralData;
  396.                         var Integral     : Float;
  397.                         var NumIntervals : integer;
  398.                         var Stack        : Ptr;
  399.                         var Finished     : boolean;
  400.                         var Error        : byte);
  401.  
  402. {-----------------------------------------------------------}
  403. {- Input: LowerLimit, UpperLimit, Tolerance, MaxIntervals  -}
  404. {- Output: Data, Integral, NumIntervals, Stack, Finished,  -}
  405. {-         Error                                           -}
  406. {-                                                         -}
  407. {- This procedure tests Tolerance and MaxIntervals for     -}
  408. {- errors (they must be greater than zero).                -}
  409. {- This procedure initializes the Data record with the     -}
  410. {- input data (LowerLimit, etc) and initializes the other  -}
  411. {- variables to 0, Nil or False.                           -}
  412. {-----------------------------------------------------------}
  413.  
  414. begin
  415.   Error := 0;
  416.   if Tolerance < 0 then
  417.     Error := 1;
  418.   if MaxIntervals < 0 then
  419.     Error := 2;
  420.   if Error = 0 then
  421.   begin
  422.     with Data do
  423.     begin
  424.       LB := LowerLimit;
  425.       VLB := UserFunction(LowerLimit, FuncPtr);
  426.       UB := UpperLimit;
  427.       VUB := UserFunction(UpperLimit, FuncPtr);
  428.       Est := Simpson(LB, VLB, UB, VUB);
  429.     end; { with }
  430.     Integral := 0;
  431.     NumIntervals := 0;
  432.     InitializeStack(Stack);
  433.     Finished := false;
  434.   end;
  435. end; { procedure TestAndInitialize }
  436.  
  437. procedure AddIntoIntegral(NewEstimate  : Float;
  438.                       var Integral     : Float;
  439.                           NumIntervals : integer;
  440.                           MaxIntervals : integer;
  441.                       var Stack        : Ptr;
  442.                       var Data         : IntegralData;
  443.                       var Finished     : boolean;
  444.                       var Error        : byte);
  445.  
  446. {-----------------------------------------------------------}
  447. {- Input: NewEstimate, NumIntervals, MaxIntervals, Stack   -}
  448. {- Output: Integral, Stack, Data, Finished, Error          -}
  449. {-                                                         -}
  450. {- This procedure adds the integral of the subinterval     -}
  451. {- (NewEstimate) into the integral of the whole interval   -}
  452. {- (Integral).  This procedure also checks to see if the   -}
  453. {- number of intervals (NumIntervals) exceeds the maximum  -}
  454. {- number of intervals allowed (MaxIntervals).  If the     -}
  455. {- Stack is empty, the integral of the whole interval has  -}
  456. {- been approximated and Found=True, otherwise information -}
  457. {- from the next subinterval is popped off the stack.      -}
  458. {-----------------------------------------------------------}
  459.  
  460. begin
  461.   Integral := Integral + NewEstimate;
  462.   if NumIntervals = MaxIntervals then
  463.     begin  { Max number of intervals exceeded }
  464.       Finished := true;
  465.       Error := 3;
  466.     end
  467.   else     { Pop next interval off stack }
  468.            { or end if Stack = NIL       }
  469.     if Stack = nil then    { For optimization purposes }
  470.                            { no function call is made }
  471.       Finished := true
  472.     else
  473.       Pop(Data, Stack);
  474. end; { procedure AddIntoIntegral }
  475.  
  476. procedure DivideInterval(Middle       : Float;
  477.                          ValueMiddle  : Float;
  478.                          NewEstimate1 : Float;
  479.                          NewEstimate2 : Float;
  480.                      var Stack        : Ptr;
  481.                      var Data         : IntegralData);
  482.  
  483. {---------------------------------------------------------------}
  484. {- Input: Middle, ValueMiddle, NewEstimate1, NewEstimate2      -}
  485. {- Output: Stack, Data,                                        -}
  486. {-                                                             -}
  487. {- This procedure stores the information from the left half    -}
  488. {- of the interval on the stack and puts the information       -}
  489. {- from the right half of the interval into the variable Data. -}
  490. {---------------------------------------------------------------}
  491.  
  492. var
  493.   StoreData : IntegralData;
  494.  
  495. begin
  496.   StoreData.LB := Data.LB;
  497.   StoreData.VLB := Data.VLB;
  498.   StoreData.UB := Middle;
  499.   StoreData.VUB := ValueMiddle;
  500.   StoreData.Est := NewEstimate1;
  501.   Push(StoreData, Stack);
  502.   Data.LB := Middle;
  503.   Data.VLB := ValueMiddle;
  504.   Data.Est := NewEstimate2;
  505. end; { procedure DivideInterval }
  506.  
  507. begin  { procedure Adaptive_Simpson }
  508.   TestAndInitialize(LowerLimit, UpperLimit, Tolerance, MaxIntervals,
  509.                     Data, Integral, NumIntervals, Stack, Finished, Error);
  510.   if Error = 0 then
  511.   repeat
  512.     with Data do
  513.     begin
  514.       NumIntervals := Succ(NumIntervals);
  515.       Spacing := (UB - LB) / 2;
  516.       Middle := LB + Spacing;
  517.       ValueMiddle := UserFunction(Middle, FuncPtr);
  518.       Estimate := Est;
  519.       { NewEstimate1 is the integral of the left half of the interval }
  520.       NewEstimate1 := Simpson(LB, VLB, Middle, ValueMiddle);
  521.       { NewEstimate2 is the integral of the right half of the interval }
  522.       NewEstimate2 := Simpson(Middle, ValueMiddle, UB, VUB);
  523.       NewEstimate := NewEstimate1 + NewEstimate2;
  524.     end;
  525.  
  526.     if (ABS(Estimate - NewEstimate) <= ABS(Tolerance * NewEstimate)) or
  527.        (NumIntervals = MaxIntervals) then { Prevents infinite loops }
  528.  
  529.       { The integral of this subinterval has      }
  530.       { been approximated to the proper tolerance }
  531.       AddIntoIntegral(NewEstimate, Integral, NumIntervals,
  532.                       MaxIntervals, Stack, Data, Finished, Error)
  533.  
  534.     else  { Divide the interval in 2 and }
  535.           { compute the integral in each }
  536.       DivideInterval(Middle, ValueMiddle, NewEstimate1, NewEstimate2,
  537.                      Stack, Data);
  538.   until Finished;
  539. end; { procedure Adaptive_Simpson }
  540.  
  541. procedure Adaptive_Gauss_Quadrature{(LowerLimit  : Float;
  542.                                     UpperLimit   : Float;
  543.                                     Tolerance    : Float;
  544.                                     MaxIntervals : integer;
  545.                                 var Integral     : Float;
  546.                                 var NumIntervals : integer;
  547.                                 var Error        : byte;
  548.                                     FuncPtr      : Pointer)};
  549.  
  550. type
  551.   IntegralData = record
  552.                    LB, VLB, UB, VUB, Est : Float;
  553.                  end;
  554.  
  555.   Ptr = ^StackItem;
  556.  
  557.   StackItem = record
  558.                 Info : IntegralData;
  559.                 Next : Ptr;
  560.               end;
  561.  
  562. var
  563.   Spacing, Estimate, NewEstimate1, NewEstimate2, NewEstimate,
  564.   Middle, ValueMiddle : Float;
  565.   Data : IntegralData;
  566.   Stack : Ptr;
  567.   Finished : boolean;
  568.  
  569. function Gaussian_Quadrature(LowerLimit : Float;
  570.                              UpperLimit : Float) : Float;
  571.  
  572. {------------------------------------------------------------}
  573. {- Input: LowerLimit, UpperLimit                            -}
  574. {-                                                          -}
  575. {- This function returns the integral of TNTargetF(X)       -}
  576. {- from LowerLimit to UpperLimit using Gaussian quadrature. -}
  577. {------------------------------------------------------------}
  578.  
  579. const
  580.   NumLegendreTerms = 16;
  581.  
  582. type
  583.   TNConstants = array[1..NumLegendreTerms] of record
  584.                                                 Root, Weight : Float;
  585.                                               end;
  586.  
  587. const
  588.  
  589. {------------------------------------------------------------------}
  590. {- These numbers are the zeros and weight factors of the          -}
  591. {- NumLegendreTermsth order Legendre Polynomial.  The numbers are -}
  592. {- taken from Abramowitz, Milton and Stegum, Irene. "Handbook of  -}
  593. {- Mathematical Functions with Formulas, Graphs and Mathematical  -}
  594. {- Tables." Washington DC: National Bureau of Standards, Applied  -}
  595. {- Mathematics Series, 55. 1972.                                  -}
  596. {-                                                                -}
  597. {- These numbers satisfy the following:                           -}
  598. {-                                                                -}
  599. {-     Integral from -1 to 1 of f(X) dx                           -}
  600. {-                  equals                                        -}
  601. {-     Sum from I=1 to NumLegendreTerms of                        -}
  602. {-                Legendre[I].Weight - f(Legendre[I].Root)        -}
  603. {------------------------------------------------------------------}
  604.  
  605.       Legendre : TNConstants = (
  606.         { Legendre[1] } (Root : 0.095012509837637440185;
  607.           Weight : 0.189450610455068496285),
  608.         { Legendre[2] } (Root : 0.281603550779258913230;
  609.           Weight : 0.182603415044923588867),
  610.         { Legendre[3] } (Root : 0.458016777657227386342;
  611.           Weight : 0.169156519395002538189),
  612.         { Legendre[4] } (Root : 0.617876244402643748447;
  613.           Weight : 0.149595988816576732081),
  614.         { Legendre[5] } (Root : 0.755404408355003033895;
  615.           Weight : 0.124628971255533872052),
  616.         { Legendre[6] } (Root : 0.865631202387831743880;
  617.           Weight : 0.095158511682492784810),
  618.         { Legendre[7] } (Root : 0.944575023073232576078;
  619.           Weight : 0.062253523938647892863),
  620.         { Legendre[8] } (Root : 0.989400934991649932596;
  621.           Weight : 0.027152459411754094852),
  622.         { Legendre[9] } (Root : -0.095012509837637440185;
  623.           Weight : 0.189450610455068496285),
  624.         { Legendre[10] } (Root : -0.281603550779258913230;
  625.           Weight : 0.182603415044923588867),
  626.         { Legendre[11] } (Root : -0.458016777657227386342;
  627.           Weight : 0.169156519395002538189),
  628.         { Legendre[12] } (Root : -0.617876244402643748447;
  629.           Weight : 0.149595988816576732081),
  630.         { Legendre[13] } (Root : -0.755404408355003033895;
  631.           Weight : 0.124628971255533872052),
  632.         { Legendre[14] } (Root : -0.865631202387831743880;
  633.           Weight : 0.095158511682492784810),
  634.         { Legendre[15] } (Root : -0.944575023073232576078;
  635.           Weight : 0.062253523938647892863),
  636.         { Legendre[16] } (Root : -0.989400934991649932596;
  637.           Weight : 0.027152459411754094852));
  638.  
  639. var
  640.   Slope, Intercept, TransRoot, Sum : Float;
  641.   Term : integer;
  642.  
  643. begin { function Gaussian_Quadrature }
  644.   Slope := (UpperLimit - LowerLimit) / 2;
  645.   Intercept := (UpperLimit + LowerLimit) / 2;
  646.   Sum := 0;
  647.   for Term := 1 to NumLegendreTerms do
  648.   with Legendre[Term] do
  649.   begin
  650.     TransRoot := Slope * Root + Intercept;
  651.     Sum := Sum + Weight * UserFunction(TransRoot, FuncPtr);
  652.   end;
  653.   Gaussian_Quadrature := Slope * Sum;
  654. end; { function Gaussian_Quadrature }
  655.  
  656. procedure InitializeStack(var Stack : Ptr);
  657. begin
  658.   Stack := nil;
  659. end; { procedure InitializeStack }
  660.  
  661. procedure Push(Data  : IntegralData;
  662.            var Stack : Ptr);
  663.  
  664. {------------------------------}
  665. {- Input: Data, Stack         -}
  666. {- Output: Stack              -}
  667. {-                            -}
  668. {- Push Data onto the Stack   -}
  669. {------------------------------}
  670.  
  671. var
  672.   NewNode : Ptr;
  673.  
  674. begin
  675.   New(NewNode);
  676.   NewNode^.Info := Data;
  677.   NewNode^.Next := Stack;
  678.   Stack := NewNode;
  679. end; { procedure Push }
  680.  
  681. procedure Pop(var Data  : IntegralData;
  682.               var Stack : Ptr);
  683.  
  684. {------------------------------}
  685. {- Input: Stack               -}
  686. {- Output: Data, Stack        -}
  687. {-                            -}
  688. {- Pop Data from the Stack    -}
  689. {------------------------------}
  690.  
  691. var
  692.   OldNode : Ptr;
  693.  
  694. begin
  695.   OldNode := Stack;
  696.   Stack := OldNode^.Next;
  697.   Data := OldNode^.Info;
  698.   Dispose(OldNode);
  699. end; { procedure Pop }
  700.  
  701. procedure TestAndInitialize(LowerLimit   : Float;
  702.                             UpperLimit   : Float;
  703.                             Tolerance    : Float;
  704.                             MaxIntervals : integer;
  705.                         var Data         : IntegralData;
  706.                         var Integral     : Float;
  707.                         var NumIntervals : integer;
  708.                         var Stack        : Ptr;
  709.                         var Finished     : boolean;
  710.                         var Error        : byte);
  711.  
  712. {-----------------------------------------------------------}
  713. {- Input: LowerLimit, UpperLimit, Tolerance, MaxIntervals  -}
  714. {- Output: Data, Integral, NumIntervals, Stack, Finished,  -}
  715. {-         Error                                           -}
  716. {-                                                         -}
  717. {- This procedure tests Tolerance and MaxIntervals for     -}
  718. {- errors (they must be greater than zero).                -}
  719. {- This procedure initializes the Data record with the     -}
  720. {- input data (LowerLimit, etc) and initializes the other  -}
  721. {- variables to 0, Nil or False.                           -}
  722. {-----------------------------------------------------------}
  723.  
  724. begin
  725.   Error := 0;
  726.   if Tolerance < 0 then
  727.     Error := 1;
  728.   if MaxIntervals < 0 then
  729.     Error := 2;
  730.   if Error = 0 then
  731.   begin
  732.     with Data do
  733.     begin
  734.       LB := LowerLimit;
  735.       VLB := UserFunction(LowerLimit, FuncPtr);
  736.       UB := UpperLimit;
  737.       VUB := UserFunction(UpperLimit, FuncPtr);
  738.       Est := Gaussian_Quadrature(LB, UB);
  739.     end; { with }
  740.     Integral := 0;
  741.     NumIntervals := 0;
  742.     InitializeStack(Stack);
  743.     Finished := false;
  744.   end;
  745. end; { procedure TestAndInitialize }
  746.  
  747. procedure AddIntoIntegral(NewEstimate  : Float;
  748.                       var Integral     : Float;
  749.                           NumIntervals : integer;
  750.                           MaxIntervals : integer;
  751.                       var Stack        : Ptr;
  752.                       var Data         : IntegralData;
  753.                       var Finished     : boolean;
  754.                       var Error        : byte);
  755.  
  756. {-----------------------------------------------------------}
  757. {- Input: NewEstimate, NumIntervals, MaxIntervals, Stack   -}
  758. {- Output: Integral, Stack, Data, Finished, Error          -}
  759. {-                                                         -}
  760. {- This procedure adds the integral of the subinterval     -}
  761. {- (NewEstimate) into the integral of the whole interval   -}
  762. {- (Integral).  This procedure also checks to see if the   -}
  763. {- number of intervals (NumIntervals) exceeds the maximum  -}
  764. {- number of intervals allowed (MaxIntervals).  If the     -}
  765. {- Stack is empty, the integral of the whole interval has  -}
  766. {- been approximated and Found=True, otherwise information -}
  767. {- from the next subinterval is popped off the stack.      -}
  768. {-----------------------------------------------------------}
  769.  
  770. begin
  771.   Integral := Integral + NewEstimate;
  772.   if NumIntervals = MaxIntervals then
  773.     begin    { Max number of intervals exceeded }
  774.       Finished := true;
  775.       Error := 3;
  776.     end
  777.   else       { Pop next interval off stack }
  778.              { or end if Stack = NIL       }
  779.     if Stack = nil then     { for optimization purposes }
  780.                             { no function call is made }
  781.       Finished := true
  782.     else
  783.       Pop(Data, Stack);
  784. end; { procedure AddIntoIntegral }
  785.  
  786. procedure DivideInterval(Middle       : Float;
  787.                          ValueMiddle  : Float;
  788.                          NewEstimate1 : Float;
  789.                          NewEstimate2 : Float;
  790.                      var Stack        : Ptr;
  791.                      var Data         : IntegralData);
  792.  
  793. {---------------------------------------------------------------}
  794. {- Input: Middle, ValueMiddle, NewEstimate1, NewEstimate2      -}
  795. {- Output: Stack, Data,                                        -}
  796. {-                                                             -}
  797. {- This procedure stores the information from the left half    -}
  798. {- of the interval on the stack and puts the information       -}
  799. {- from the right half of the interval into the variable Data. -}
  800. {---------------------------------------------------------------}
  801.  
  802. var
  803.   StoreData : IntegralData;
  804.  
  805. begin
  806.   StoreData.LB := Data.LB;
  807.   StoreData.VLB := Data.VLB;
  808.   StoreData.UB := Middle;
  809.   StoreData.VUB := ValueMiddle;
  810.   StoreData.Est := NewEstimate1;
  811.   Push(StoreData, Stack);
  812.   Data.LB := Middle;
  813.   Data.VLB := ValueMiddle;
  814.   Data.Est := NewEstimate2;
  815. end; { procedure DivideInterval }
  816.  
  817. begin  { procedure Adaptive_Gauss_Quadrature }
  818.   TestAndInitialize(LowerLimit, UpperLimit, Tolerance, MaxIntervals,
  819.                     Data, Integral, NumIntervals, Stack, Finished, Error);
  820.  
  821.   if Error = 0 then
  822.   repeat
  823.     with Data do
  824.     begin
  825.       NumIntervals := Succ(NumIntervals);
  826.       Spacing := (UB - LB) / 2;
  827.       Middle := LB + Spacing;
  828.       ValueMiddle := UserFunction(Middle, FuncPtr);
  829.       Estimate := Est;
  830.       { NewEstimate1 is the integral of the left half of the interval }
  831.       NewEstimate1 := Gaussian_Quadrature(LB, Middle);
  832.       { NewEstimate2 is the integral of the right half of the interval }
  833.       NewEstimate2 := Gaussian_Quadrature(Middle, UB);
  834.       NewEstimate := NewEstimate1 + NewEstimate2;
  835.     end; { with }
  836.  
  837.     if (ABS(Estimate - NewEstimate) <= ABS(Tolerance * NewEstimate)) or
  838.        (NumIntervals = MaxIntervals) then { Prevents infinite loops }
  839.  
  840.       { The integral of this subinterval has      }
  841.       { been approximated to the proper tolerance }
  842.       AddIntoIntegral(NewEstimate, Integral, NumIntervals,
  843.                       MaxIntervals, Stack, Data, Finished, Error)
  844.  
  845.     else  { Divide the interval in 2 and }
  846.           { compute the integral in each }
  847.       DivideInterval(Middle, ValueMiddle, NewEstimate1, NewEstimate2,
  848.                      Stack, Data);
  849.   until Finished;
  850. end; { procedure Adaptive_Gauss_Quadrature }
  851.  
  852. procedure Romberg{(LowerLimit : Float;
  853.                   UpperLimit  : Float;
  854.                   Tolerance   : Float;
  855.                   MaxIter     : integer;
  856.               var Integral    : Float;
  857.               var Iter        : integer;
  858.               var Error       : byte;
  859.                   FuncPtr     : Pointer)};
  860.  
  861. var
  862.   Spacing : Float;                 { Spacing between points }
  863.   NewEstimate,
  864.   OldEstimate : TNvector;         { Iteration variables }
  865.   TwoToTheIterMinus2 : integer;
  866.  
  867. procedure TestAndInitialize(LowerLimit         : Float;
  868.                             UpperLimit         : Float;
  869.                             Tolerance          : Float;
  870.                             MaxIter            : integer;
  871.                         var Iter               : integer;
  872.                         var Spacing            : Float;
  873.                         var OldEstimate        : TNvector;
  874.                         var TwoToTheIterMinus2 : integer;
  875.                         var Error              : byte);
  876.  
  877. {------------------------------------------------------------------}
  878. {- Input: LowerLimit, UpperLimit, Tolerance, MaxIter              -}
  879. {- Output: Iter, Spacing, OldEstimate, TwoToTheIterMinus2, Error  -}
  880. {-                                                                -}
  881. {- This procedure tests Tolerance and MaxIter for errors (they    -}
  882. {- must be greater than zero) and initializes the above           -}
  883. {- variables.                                                     -}
  884. {------------------------------------------------------------------}
  885.  
  886. begin
  887.   Error := 0;
  888.   if Tolerance <= 0 then
  889.     Error := 1;
  890.   if MaxIter <= 0 then
  891.     Error := 2;
  892.   if Error = 0 then
  893.   begin
  894.     Spacing := UpperLimit - LowerLimit;
  895.     OldEstimate[1] := Spacing *
  896.                       (UserFunction(LowerLimit, FuncPtr) +
  897.                        UserFunction(UpperLimit, FuncPtr)) / 2;
  898.     Iter := 1;
  899.     TwoToTheIterMinus2 := 1;
  900.   end;
  901. end; { procedure TestAndInitialize }
  902.  
  903. procedure Trapezoid(TwoToTheIterMinus2 : integer;
  904.                     LowerLimit         : Float;
  905.                     Spacing            : Float;
  906.                     OldEstimate        : Float;
  907.                 var NewEstimate        : Float);
  908.  
  909. {--------------------------------------------------------}
  910. {- Input: TwoToTheIterMinus2, LowerLimit, Spacing,      -}
  911. {-        OldEstimate                                   -}
  912. {- Output: NewEstimate                                  -}
  913. {-                                                      -}
  914. {- This procedure uses the trapezoid rule to            -}
  915. {- improve the integral approximation (OldEstimate)     -}
  916. {- on the interval [LowerLimit, LowerLimit + Spacing].  -}
  917. {- The results are returned in the variable NewEstimate -}
  918. {--------------------------------------------------------}
  919.  
  920. var
  921.   Sum : Float;
  922.   Dummy : integer;
  923.  
  924. begin
  925.   Sum := 0;
  926.   for Dummy := 1 to TwoToTheIterMinus2 do
  927.     Sum := Sum + UserFunction(LowerLimit + (Dummy - 0.5) * Spacing, FuncPtr);
  928.   NewEstimate := 0.5 * (OldEstimate + Spacing * Sum);
  929. end; { procedure Trapezoid }
  930.  
  931. procedure Extrapolate(Iter        : integer;
  932.                       OldEstimate : TNvector;
  933.                   var NewEstimate : TNvector);
  934.  
  935. {--------------------------------------------------------}
  936. {- Input: Iter, OldEstimate                             -}
  937. {- Output: NewEstimate                                  -}
  938. {-                                                      -}
  939. {- This procedure uses Richardson extrapolation         -}
  940. {- to improve the current approximation to the integral -}
  941. {- (OldEstimate).  The result is returned in the        -}
  942. {- variable NewEstimate                                 -}
  943. {--------------------------------------------------------}
  944.  
  945. var
  946.   Extrap : integer;
  947.   FourToTheExtrapMinus1 : Float;
  948.  
  949. begin
  950.   FourToTheExtrapMinus1 := 1;
  951.   for Extrap := 2 to Iter do
  952.   begin
  953.     FourToTheExtrapMinus1 := FourToTheExtrapMinus1 * 4;
  954.     NewEstimate[Extrap] :=
  955.                         (FourToTheExtrapMinus1 * NewEstimate[Extrap - 1] -
  956.                         OldEstimate[Extrap - 1]) / (FourToTheExtrapMinus1 - 1);
  957.   end;
  958. end; { procedure Extrapolate }
  959.  
  960. begin { procedure Romberg }
  961.   TestAndInitialize(LowerLimit, UpperLimit, Tolerance, MaxIter, Iter,
  962.                     Spacing, OldEstimate, TwoToTheIterMinus2, Error);
  963.   if Error = 0 then
  964.   begin
  965.     repeat
  966.       Iter := Succ(Iter);
  967.       Trapezoid(TwoToTheIterMinus2, LowerLimit, Spacing, OldEstimate[1],
  968.                 NewEstimate[1]);
  969.       TwoToTheIterMinus2 := TwoToTheIterMinus2 * 2;
  970.       Extrapolate(Iter, OldEstimate, NewEstimate);
  971.       Spacing := Spacing / 2;
  972.       OldEstimate := NewEstimate;
  973.     until { The fractional difference between }
  974.           { iterations is within Tolerance    }
  975.           (ABS(NewEstimate[Iter - 1] - NewEstimate[Iter]) <=
  976.            ABS(Tolerance * NewEstimate[Iter])) or (Iter >= MaxIter);
  977.  
  978.     if Iter >= MaxIter then
  979.       Error := 3;
  980.     Integral := NewEstimate[Iter];
  981.   end;
  982. end; { procedure Romberg }
  983.  
  984. end. { Integrat }
  985.