home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / pascal / library / dos / sampler / 01 / num_tool / raphson.inc < prev    next >
Encoding:
Text File  |  1987-10-22  |  8.9 KB  |  188 lines
  1. PROCEDURE Newton_Raphson(Guess : Real;
  2.                          Tol : Real;
  3.                          MaxIter : Integer;
  4.                          VAR Root : Real;
  5.                          VAR Value : Real;
  6.                          VAR Deriv : Real;
  7.                          VAR Iter : Integer;
  8.                          VAR Error : Byte);
  9.  
  10. {------------------------------------------------------------------}
  11. {-                                                                -}
  12. {- Turbo Pascal Numerical Methods Toolbox                         -}
  13. {- (C) Copyright 1986 Borland International.                      -}
  14. {-                                                                -}
  15. {-         Input: Guess, Tol, MaxIter                             -}
  16. {-         Output: Root, Value, Deriv, Iter, Error                -}
  17. {-                                                                -}
  18. {- Purpose: This unit provides a procedure for finding a single   -}
  19. {-          real root of a user specified function with a known   -}
  20. {-          continuous first derivative, given a user             -}
  21. {-          specified initial guess.  The procedure implements    -}
  22. {-          Newton-Raphson's algorithm for finding a single       -}
  23. {-          zero of a function.                                   -}
  24. {-          The user must specify the desired tolerance           -}
  25. {-          in the answer.                                        -}
  26. {-                                                                -}
  27. {- Global Variables:                                              -}
  28. {-    Guess   : real;    user's estimate of root                  -}
  29. {-            Tol     : real;    tolerance in answer              -}
  30. {-            MaxIter : integer; maximum number of iterations     -}
  31. {-            Root    : real;    real part of calculated roots    -}
  32. {-            Value   : real;    value of the polynomial at root  -}
  33. {-            Deriv   : real;    value of the derivative at root  -}
  34. {-            Iter    : real;    number of iterations it took     -}
  35. {-                               to find root                     -}
  36. {-            Error   : byte;    flags if something went wrong    -}
  37. {-                                                                -}
  38. {-    Errors: 1: Iter >= MaxIter                                  -}
  39. {-            2: The slope was zero at some point                 -}
  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, you are }
  49.   { not running TURBO-87.                      }
  50.   { TNNearlyZero = 1E-015 if using the 8087    }
  51.   {                      math co-processor.    }
  52.   { TNNearlyZero = 1E-07 if not using the 8087 }
  53.   {                         math co-processor. }
  54.  
  55. VAR
  56.   Found : Boolean;          { Flags that a root has been Found }
  57.   OldX, OldY, OldDeriv,
  58.   NewX, NewY, NewDeriv : Real; { Iteration variables }
  59.  
  60.   PROCEDURE CheckSlope(Slope : Real;
  61.                        VAR Error : Byte);
  62.  
  63.     {---------------------------------------------------}
  64.     {- Input:  Slope                                   -}
  65.     {- Output: Error                                   -}
  66.     {-                                                 -}
  67.     {- This procedure checks the slope to see if it is -}
  68.     {- zero.  The Newton Raphson algorithm may not be  -}
  69.     {- applied at a point where the slope is zero.     -}
  70.     {---------------------------------------------------}
  71.  
  72.   BEGIN
  73.     IF Abs(Slope) <= TNNearlyZero THEN
  74.       Error := 2;           { Slope is zero }
  75.   END;                      { procedure CheckSlope }
  76.  
  77.   PROCEDURE Initial(Guess : Real;
  78.                     Tol : Real;
  79.                     MaxIter : Integer;
  80.                     VAR OldX : Real;
  81.                     VAR OldY : Real;
  82.                     VAR OldDeriv : Real;
  83.                     VAR Found : Boolean;
  84.                     VAR Iter : Integer;
  85.                     VAR Error : Byte);
  86.  
  87.  
  88.     {-------------------------------------------------------------}
  89.     {- Input:  Guess, Tol, MaxIter                               -}
  90.     {- Output: OldX, OldY, OldDeriv, Found, Iter, Error          -}
  91.     {-                                                           -}
  92.     {- This procedure sets the initial values of the above       -}
  93.     {- variables. If OldY is zero, then a root has been          -}
  94.     {- Found and Found = TRUE.  This procedure also checks       -}
  95.     {- the tolerance (Tol) and the maximum number of iterations  -}
  96.     {- (MaxIter) for errors.                                     -}
  97.     {-------------------------------------------------------------}
  98.  
  99.   BEGIN
  100.     Found := False;
  101.     Iter := 0;
  102.     Error := 0;
  103.     OldX := Guess;
  104.     OldY := TNTargetF(OldX);
  105.     OldDeriv := TNDerivF(OldX);
  106.     IF Abs(OldY) <= TNNearlyZero THEN
  107.       Found := True
  108.     ELSE
  109.       CheckSlope(OldDeriv, Error);
  110.     IF Tol <= 0 THEN
  111.       Error := 3;
  112.     IF MaxIter < 0 THEN
  113.       Error := 4;
  114.   END;                      { procedure Initial }
  115.  
  116.   FUNCTION TestForRoot(X, OldX, Y, Tol : Real) : Boolean;
  117.  
  118. {-----------------------------------------------------------------}
  119. {- These are the stopping criteria.  Four different ones are     -}
  120. {- provided.  If you wish to change the active criteria, simply  -}
  121. {- comment off the current criteria (including the preceding OR) -}
  122. {- and remove the comment brackets from the criteria (including  -}
  123. {- the following OR) you wish to be active.                      -}
  124. {-----------------------------------------------------------------}
  125.  
  126.   BEGIN
  127.  
  128.     TestForRoot :=                    {---------------------------}
  129.       (ABS(Y) <= TNNearlyZero)        {- Y = 0                   -}
  130.                                       {-                         -}
  131.              or                       {-                         -}
  132.                                       {-                         -}
  133.       (ABS(X - OldX) < ABS(OldX*Tol)) {- Relative change in X    -}
  134.                                       {-                         -}
  135. (*           or                    *) {-                         -}
  136. (*                                 *) {-                         -}
  137. (*    (ABS(OldX - X) < Tol)        *) {- Absolute change in X    -}
  138. (*                                 *) {-                         -}
  139. (*           or                    *) {-                         -}
  140. (*                                 *) {-                         -}
  141. (*    (ABS(Y) <= Tol)              *) {- Absolute change in Y    -}
  142.                                       {---------------------------}
  143.  
  144. {-------------------------------------------------------------------}
  145. {- The first criteria simply checks to see if the value of the     -}
  146. {- function is zero.  You should probably always keep this         -}
  147. {- criteria active.                                                -}
  148. {-                                                                 -}
  149. {- The second criteria checks relative error in X. This criteria   -}
  150. {- evaluates the fractional change in X between interations.       -}
  151. {- Note that X has been multiplied throught the inequality to      -}
  152. {- avoid divide by zero errors.                                    -}
  153. {-                                                                 -}
  154. {- The third criteria checks the absolute difference in X          -}
  155. {- between iterations.                                             -}       
  156. {-                                                                 -}  
  157. {- The fourth criteria checks the absolute difference between      -}
  158. {- the value of the function and zero.                             -}
  159. {-                                                                 -}
  160. {-------------------------------------------------------------------}
  161.  
  162.   END;                      { procedure TestForRoot }
  163.  
  164.  
  165. BEGIN                       { procedure Newton_Raphson }
  166.   Initial(Guess, Tol, MaxIter, OldX, OldY, OldDeriv,
  167.   Found, Iter, Error);
  168.  
  169.   WHILE NOT(Found) AND (Error = 0) AND (Iter < MaxIter) DO
  170.     BEGIN
  171.       Iter := Succ(Iter);
  172.       NewX := OldX-OldY/OldDeriv;
  173.       NewY := TNTargetF(NewX);
  174.       NewDeriv := TNDerivF(NewX);
  175.       Found := TestForRoot(NewX, OldX, NewY, Tol);
  176.       OldX := NewX;
  177.       OldY := NewY;
  178.       OldDeriv := NewDeriv;
  179.       IF NOT(Found) THEN
  180.         CheckSlope(OldDeriv, Error);
  181.     END;
  182.   Root := OldX;
  183.   Value := OldY;
  184.   Deriv := OldDeriv;
  185.   IF NOT(Found) AND (Error = 0) AND (Iter >= MaxIter) THEN
  186.     Error := 1;
  187. END;                        { procedure Newton_Raphson }
  188.