home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l292 / 1.ddi / ROOTS.FOR < prev    next >
Encoding:
Text File  |  1989-10-10  |  5.6 KB  |  221 lines
  1.  
  2.       REAL FUNCTION AFunction (X)
  3.       REAL X
  4.         AFunction = (2 * X + 3) * (X - 3)
  5.       RETURN
  6.       END !FUNCTION
  7.  
  8.       REAL FUNCTION AFunctionPrime (X)
  9.       REAL X
  10.         AFunctionPrime = 4 * X - 3
  11.       RETURN
  12.       END  !FUNCTION
  13.  
  14.  
  15.  
  16.  
  17.       FUNCTION BisectionRoots (x1, x2, Tolerance, maxIterations,
  18.      +        valueAtRoot, errr)
  19.       REAL x1, x2, Tolerance, valueAtRoot, FPP, nearzero
  20.       INTEGER maxIterations, errr
  21.  
  22.       PARAMETER (FPP = 1E-15)
  23.       PARAMETER (nearzero = 1E-20 )
  24.       REAL result, FA,  FC, XA, XB, XC
  25.       INTEGER i
  26.       LOGICAL done
  27.  
  28.  
  29.       errr = 0
  30.       i = 0
  31.       done = .FALSE.
  32.       XA = x1
  33.       XB = x2
  34.       FA = AFunction(XA)
  35.       FB = AFunction(XB)
  36.       IF ( RootBracketed(FA, FB) .EQ. 0 ) THEN
  37.         errr = 1
  38.       ELSE
  39.         DO WHILE (.NOT. done .AND. i .NE. maxIterations)
  40.           XC = (XA + XB) / 2.0
  41.           FC = AFunction(XC)
  42.           IF (RootBracketed(FA, FC) .EQ. 1) THEN
  43.             XB = XC
  44.           ELSE
  45.             XA = XC
  46.             FA = FC
  47.           END IF
  48.           IF (ABS(XB - XA) .LE. Tolerance) THEN
  49.             done = .TRUE.
  50.             result = (XA + XB) / 2.0
  51.             valueAtRoot = AFunction(result)
  52.           END IF
  53.           i = i + 1
  54.         END DO
  55.         IF (i .EQ. maxIterations) errr = 2
  56.       END IF
  57.       BisectionRoots = result
  58.       END !FUNCTION
  59.  
  60.  
  61.  
  62.       REAL FUNCTION BrentRoots (x1, x2, Tolerance, maxIterations,
  63.      +    valueAtRoot, errr)
  64.       REAL x1, x2, Tolerance, valueAtRoot, FPP, nearzero
  65.       INTEGER maxIterations, errr
  66.  
  67.       PARAMETER (FPP = 1E-15)
  68.       PARAMETER (nearzero = 1E-20)
  69.       REAL result, AA, BB, CC, DD, EE, FA, FB
  70.       REAL  Tol1, PP, QQ, RR, SS, XM
  71.       INTEGER i
  72.       LOGICAL done
  73.       i = 0
  74.       done = .FALSE.
  75.       errr = 0
  76.       AA = x1
  77.       BB = x2
  78.       FA = AFunction(AA)
  79.       FB = AFunction(BB)
  80.       IF ( RootBracketed(FA, FB) .EQ. 0) THEN
  81.          errr = 1
  82.       ELSE
  83.         FC = FB
  84.         DO WHILE (.NOT. done .AND. i .NE. maxIterations)
  85.           IF  (RootBracketed(FC, FB) .EQ. 0) THEN
  86.             CC = AA
  87.             FC = FA
  88.             DD = BB - AA
  89.             EE = DD
  90.           END IF
  91.           IF (ABS(FC) .LT. ABS(FB)) THEN
  92.             AA = BB
  93.             BB = CC
  94.             CC = AA
  95.             FA = FB
  96.             FB = FC
  97.             FC = FA
  98.           END IF
  99.           Tol1 = 2.0 * FPP * ABS(BB) + 0.5 * Tolerance
  100.           XM = 0.5 * (CC - BB)
  101.           IF ((ABS(XM) .LE. Tol1) .OR. (ABS(FA) .LT. nearzero)) THEN
  102.             result = BB
  103.             done = .TRUE.
  104.             valueAtRoot = AFunction(result)
  105.           ELSE
  106.             IF ((ABS(EE) .GE. Tol1) .AND. (ABS(FA) .GT. ABS(FB))) THEN
  107.               SS = FB / FA
  108.               IF (ABS(AA - CC) .LT. nearzero) THEN
  109.                 PP = 2.0 * XM * SS
  110.                 QQ = 1.0 - SS
  111.               ELSE
  112.                 QQ = FA / FC
  113.                 RR = FB / FC
  114.                 PP = SS * (2.0 * XM * QQ *
  115.      +               (QQ - RR) - (BB - AA) * (RR - 1.0))
  116.                 QQ = (QQ - 1.0) * (RR - 1.0) * (SS - 1.0)
  117.               END IF
  118.               IF (PP .GT. nearzero)  QQ = -QQ
  119.               PP = ABS(PP)
  120.               IF (2.0 * PP .LT.
  121.      +           Minimum(3.0*XM*QQ - ABS(Tol1*QQ), ABS(EE*QQ))) THEN
  122.                 EE = DD
  123.                 DD = PP / QQ
  124.               ELSE
  125.                 DD = XM
  126.                 EE = DD
  127.               END IF
  128.             ELSE
  129.               DD = XM
  130.               EE = DD
  131.             END IF
  132.             AA = BB
  133.             FA = FB
  134.             IF (ABS(DD) .GT. Tol1) THEN
  135.               BB = BB + DD
  136.             ELSE
  137.               IF (XM .GT. 0) THEN
  138.                 BB = BB + ABS(Tol1)
  139.               ELSE
  140.                 BB = BB - ABS(Tol1)
  141.               END IF
  142.             END IF
  143.             FB = AFunction(BB)
  144.           END IF
  145.           i = i + 1
  146.         END DO
  147.         IF (i .EQ. maxIterations)  errr = 2
  148.       END IF
  149.       BrentRoots = result
  150.       END !FUNCTION
  151.  
  152.  
  153.  
  154.       FUNCTION Minimum (x1, x2)
  155.       REAL x1, x2
  156.       REAL result
  157.  
  158.       IF (x1 .LT. x2) THEN
  159.         result = x1
  160.       ELSE
  161.         result = x2
  162.       END IF
  163.       Minimum = result
  164.       END !FUNCTION
  165.  
  166.  
  167.  
  168.       FUNCTION NewtonRoots (x1, Tolerance, maxIterations,
  169.      +     valueAtRoot, errr)
  170.       REAL FPP, nearzero
  171.       PARAMETER ( FPP = 1E-15)
  172.       PARAMETER ( nearzero = 1E-20)
  173.       REAL result, lastX, FX, FPX, X
  174.       LOGICAL done
  175.       INTEGER maxIterations, errr
  176.  
  177.       i = 0
  178.       errr = 0
  179.       done = .FALSE.
  180.       X = x1
  181.       FX = AFunction(X)
  182.       DO WHILE (.NOT. done .AND. i .LE. maxIterations .AND. errr .EQ. 0)
  183.  
  184.        IF (ABS(X) .LT. nearzero) THEN
  185.          errr = 1
  186.        ELSE
  187.          FPX = AFunctionPrime(X)
  188.          IF (ABS(FPX) .LT. nearzero) THEN
  189.            errr = 2
  190.          ELSE
  191.            lastX = X
  192.            X = X - FX / FPX
  193.            FX = AFunction(X)
  194.            IF (ABS(X - lastX) .LE. Tolerance) THEN
  195.              done = .TRUE.
  196.              result = X
  197.              valueAtRoot = AFunction(X)
  198.            END IF
  199.            i = i + 1
  200.          END IF
  201.         END IF
  202.       END DO !! UNTIL done .OR. (i .GT. maxIterations) .OR. (errr .NE. 0)
  203.       IF( i .EQ. maxIterations)  errr = 3
  204.       NewtonRoots = result
  205.  
  206.       END ! FUNCTION
  207.  
  208.  
  209.       FUNCTION RootBracketed (x1, x2)
  210.       REAL x1, x2
  211.       REAL result
  212.       IF ( (x1 .GT. 0.0 .AND. x2 .GT. 0.0) .OR.
  213.      +   (x1 .LT. 0.0 .AND. x2 .LT. 0.0)) THEN
  214.         result = 0
  215.       ELSE
  216.         result =  1
  217.       END IF
  218.       RootBracketed = result
  219.       END !FUNCTION
  220.  
  221.