home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / MAGAZINE / MISC / JNFB88.ZIP / SKYDIV.ARC / BISECT.INC next >
Encoding:
Text File  |  1987-07-17  |  10.0 KB  |  224 lines
  1.  
  2. PROCEDURE Bisect(LeftEnd : Real;
  3.                  RightEnd : Real;
  4.                  Tol : Real;
  5.                  MaxIter : Integer;
  6.                  VAR Answer : Real;
  7.                  VAR fAnswer : Real;
  8.                  VAR Iter : Integer;
  9.                  VAR Error : Byte);
  10.  
  11.   {----------------------------------------------------------------}
  12.   {-                                                              -}
  13.   {-     Turbo Pascal Numerical Methods Toolbox                   -}
  14.   {-     (C) Copyright 1986 Borland International.                -}
  15.   {-                                                              -}
  16.   {- Input:  LeftEnd, RightEnd, Tol, MaxIter                      -}
  17.   {- Output: Answer, fAnswer, Iter, Error                         -}
  18.   {-                                                              -}
  19.   {- Purpose: This unit provides a procedure for finding a root   -}
  20.   {-       of a user specified function, for a user specified     -}
  21.   {-       interval, [a,b], where f(a) and f(b) are of opposite   -}
  22.   {-       signs.  The algorithm successively bisects the         -}
  23.   {-       interval, closing in on the root.  The user must       -}
  24.   {-       supply the desired tolerance to which the root should  -}
  25.   {-       be found.                                              -}
  26.   {-                                                              -}
  27.   {- Global Variables: LeftEnd  : real    left endpoint           -}
  28.   {-                   RightEnd : real    right endpoint          -}
  29.   {-                   Tol      : real    tolerance of error      -}
  30.   {-                   MaxIter  : real    max number iterations   -}
  31.   {-                   Answer   : real    root of TNTargetF       -}
  32.   {-                   fAnswer  : real    TNTargetF(Answer)       -}
  33.   {-                                      (should be close to 0)  -}
  34.   {-                   Iter     :integer  number of iterations    -}
  35.   {-                   Error    : byte    error flags             -}
  36.   {-                                                              -}
  37.   {-           Errors: 0: No error                                -}
  38.   {-                   1: maximum number of iterations exceeded   -}
  39.   {-                   2: f(a) and f(b) are not of opposite signs -}
  40.   {-                   3: Tol <= 0                                -}
  41.   {-                   4: MaxIter < 0                             -}
  42.   {-                                                              -}
  43.   {-    Version Date: 26 January 1987                             -}
  44.   {-                                                              -}
  45.   {----------------------------------------------------------------}
  46.  
  47. CONST
  48.   TNNearlyZero = 1E-015;    { If you get a syntax error here, }
  49.         
  50.   { you are not running TURBO-87.              }
  51.   { TNNearlyZero = 1E-015 if using the 8087    }
  52.   {                      math co-processor.    }
  53.   { TNNearlyZero = 1E-07 if not using the 8087 }
  54.   {                      math co-processor.    }
  55.  
  56. VAR
  57.   Found : Boolean;
  58.  
  59.   PROCEDURE TestInput(LeftEndpoint : Real;
  60.                       RightEndpoint : Real;
  61.                       Tol : Real;
  62.                       MaxIter : Integer;
  63.                       VAR Answer : Real;
  64.                       VAR fAnswer : Real;
  65.                       VAR Error : Byte;
  66.                       VAR Found : Boolean);
  67.  
  68.     {------------------------------------------------------------}
  69.     {- Input:  LeftEndpoint, RightEndpoint, Tol, MaxIter        -}
  70.     {- Output: Answer, fAnswer, Error, Found                    -}
  71.     {-                                                          -}
  72.     {- This procedure tests the input data for errors.  If      -}
  73.     {- LeftEndpoint > RightEndpoint, Tol <= 0, or MaxIter < 0,  -}
  74.     {- then an error is returned.  If one the of the endpoints  -}
  75.     {- (LeftEndpoint, RightEndpoint) is a root, then Found=TRUE -}
  76.     {- and Answer and fAnswer are returned.                     -}
  77.     {------------------------------------------------------------}
  78.  
  79.   VAR
  80.     yLeft, yRight : Real; { The values of function at endpoints. }
  81.  
  82.   BEGIN
  83.     yLeft := TNTargetF(LeftEndpoint);
  84.     yRight := TNTargetF(RightEndpoint);
  85.  
  86.     IF Abs(yLeft) <= TNNearlyZero THEN
  87.       BEGIN
  88.         Answer := LeftEndpoint;
  89.         fAnswer := TNTargetF(Answer);
  90.         Found := True;
  91.       END;
  92.  
  93.     IF Abs(yRight) <= TNNearlyZero THEN
  94.       BEGIN
  95.         Answer := RightEndpoint;
  96.         fAnswer := TNTargetF(Answer);
  97.         Found := True;
  98.       END;
  99.  
  100.     IF NOT Found THEN       { Test for errors }
  101.       BEGIN
  102.         IF yLeft*yRight > 0 THEN
  103.           Error := 2;
  104.         IF Tol <= 0 THEN
  105.           Error := 3;
  106.         IF MaxIter < 0 THEN
  107.           Error := 4;
  108.       END;
  109.   END;                      { procedure Tests }
  110.  
  111.   PROCEDURE Converge(VAR LeftEndpoint : Real;
  112.                      VAR RightEndpoint : Real;
  113.                      Tol : Real;
  114.                      VAR Found : Boolean;
  115.                      MaxIter : Integer;
  116.                      VAR Answer : Real;
  117.                      VAR fAnswer : Real;
  118.                      VAR Iter : Integer;
  119.                      VAR Error : Byte);
  120.  
  121.     {-------------------------------------------------------------}
  122.     {- Input:  LeftEndpoint, RightEndpoint, Tol, MaxIter         -}
  123.     {- Output: Found, Answer, fAnswer, Iter, Error               -}
  124.     {-                                                           -}
  125.     {- This procedure applies the bisection method to find a     -}
  126.     {- root to TNTargetF(x) on the interval [LeftEndpoint,       -}
  127.     {- RightEndpoint]. The root must be found within MaxIter     -}
  128.     {- iterations to a tolerance of Tol. If root found, then it  -}
  129.     {- is returned in Answer, the value of the function at the   -}
  130.     {- approximated root is fAnswer (should be close to          -}
  131.     {- zero), and the number of iterations is returned in Iter.  -}
  132.     {-------------------------------------------------------------}
  133.  
  134.   VAR
  135.     yLeft : Real;
  136.     MidPoint, yMidPoint : Real;
  137.  
  138.     PROCEDURE Initial(LeftEndpoint : Real;
  139.                       VAR Iter : Integer;
  140.                       VAR yLeft : Real);
  141.       { Initialize variables. }
  142.     BEGIN
  143.       Iter := 0;
  144.       yLeft := TNTargetF(LeftEndpoint);
  145.     END;                    { procedure Initial }
  146.  
  147.  
  148. FUNCTION TestForRoot(X, OldX, Y, Tol : Real) : Boolean;
  149. {-----------------------------------------------------------------}
  150. {-  These are the stopping criteria.  Four different ones are    -}
  151. {-  provided.  If you wish to change the active criteria, simply -}
  152. {-  comment off current criteria (including the appropriate OR)  -}
  153. {-  and remove the comment brackets from the criteria (including -}
  154. {-  the appropriate OR) you wish to be active.                   -}
  155. {-----------------------------------------------------------------}
  156. BEGIN
  157.   TestForRoot :=                      {------------------------}
  158.     (ABS(Y) <= TNNearlyZero)          {- Y = 0                -}
  159.                                       {-                      -}
  160.            OR                         {-                      -}
  161.                                       {-                      -}
  162.     (ABS(X - OldX) < ABS(OldX * Tol)) {- Relative change in X -}
  163.                                       {-                      -}
  164.                                       {-                      -}
  165. (*         OR                      *) {-                      -}
  166. (*                                 *) {-                      -}
  167. (*  (ABS(X - OldX) < Tol)          *) {- Absolute change in X -}
  168. (*                                 *) {-                      -}
  169. (*         OR                      *) {-                      -}
  170. (*                                 *) {-                      -}
  171. (*  (ABS(Y) <= Tol)                *) {- Absolute change in Y -}
  172.                                       {------------------------}
  173.  
  174. {------------------------------------------------------------------}
  175. {- The first criteria simply checks to see if the value of the    -}
  176. {- function is zero. You should probably always keep this         -}
  177. {- criteria  active.                                              -}
  178. {-                                                                -}
  179. {- The second criteria checks relative error in x. This criteria  -}
  180. {- evaluates the fractional change in x between interations. Note -}
  181. {- x has been multiplied through the inequality to avoid divide   -}
  182. {- by zero errors.                                                -}
  183. {-                                                                -}
  184. {- The third criteria checks the absolute difference in x         -}
  185. {- between iterations.                                            -}
  186. {-                                                                -}
  187. {- The fourth criteria checks the absolute difference between     -}
  188. {- the value of the function and zero.                            -}
  189. {------------------------------------------------------------------}
  190.  
  191. END; { function TestForRoot }
  192.  
  193.  
  194.   BEGIN                     { procedure Converge }
  195.     Initial(LeftEndpoint, Iter, yLeft);
  196.     WHILE NOT(Found) AND (Iter < MaxIter) DO
  197.       BEGIN
  198.         Iter := Succ(Iter);
  199.         MidPoint := (LeftEndpoint+RightEndpoint)/2;
  200.         yMidPoint := TNTargetF(MidPoint);
  201.         Found :=TestForRoot(MidPoint, LeftEndpoint, yMidPoint, Tol);
  202.         IF (yLeft*yMidPoint) < 0 THEN
  203.           RightEndpoint := MidPoint
  204.         ELSE
  205.           BEGIN
  206.             LeftEndpoint := MidPoint;
  207.             yLeft := yMidPoint;
  208.           END;
  209.       END;
  210.     Answer := MidPoint;
  211.     fAnswer := yMidPoint;
  212.     IF Iter >= MaxIter THEN
  213.       Error := 1;
  214.   END;                      { procedure Converge }
  215.  
  216. BEGIN                       { procedure Bisect }
  217.   Found := False;
  218.   TestInput(LeftEnd, RightEnd, Tol, MaxIter, Answer,
  219.   fAnswer, Error, Found);
  220.   IF (Error = 0) AND (Found = False) THEN { i.e. no error }
  221.     Converge(LeftEnd, RightEnd, Tol, Found, MaxIter,
  222.     Answer, fAnswer, Iter, Error);
  223. END;                        { procedure Bisect }
  224.