home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l292 / 1.ddi / INTEGRAT.FOR < prev    next >
Encoding:
Text File  |  1990-02-20  |  5.3 KB  |  189 lines
  1.  
  2.       FUNCTION AFunction (y)
  3.       REAL y
  4.         AFunction = EXP(y)
  5.       END !FUNCTION
  6.  
  7.  
  8.  
  9.       SUBROUTINE IntegrateFunction (numseg, xl, xh, ans)
  10.       INTEGER numseg
  11.       REAL xl, xh, ans
  12.       REAL strtpnt, endpnt, area1, area2, area, segwidth
  13.  
  14.       segwidth = (xh - xl) / (numseg * 1.0)
  15.  
  16.       endpnt = xh
  17.       area1 = 0.0
  18.       area2 = 0.0
  19.       area = 0.0
  20.  
  21.       IF (MOD(numseg, 2) .NE. 0) THEN
  22.       !!  SIMPSONS 3/8 rule
  23.         area1 = 3.0 / 8.0 * segwidth *
  24.      +          (AFunction(endpnt - 3 * segwidth) + 3 *
  25.      +           AFunction(endpnt - 2 * segwidth) + 3 *
  26.      +            AFunction(endpnt - segwidth) + AFunction(endpnt))
  27.         endpnt = endpnt - 3 * segwidth
  28.       ELSE
  29.         area1 = 0.0
  30.       END IF
  31.       IF (numseg .NE. 3) THEN
  32.       !!  SIMPSONS 1/3 rule
  33.         strtpnt = xl
  34.         IF  (strtpnt .GE. endpnt - segwidth) THEN
  35.           area2 = area2 + 1.0 / 3.0 * (segwidth * (
  36.      +            AFunction(strtpnt) + 4 *
  37.      +            AFunction(strtpnt + segwidth) +
  38.      +            AFunction(strtpnt + 2 * segwidth)))
  39.         strtpnt = strtpnt + 2 * segwidth
  40.         ELSE
  41.           DO WHILE (strtpnt .LT. endpnt - segwidth)
  42.             area2 = area2 + 1.0 / 3.0 * (segwidth * (
  43.      +              AFunction(strtpnt) + 4 *
  44.      +              AFunction(strtpnt + segwidth) +
  45.      +              AFunction(strtpnt + 2 * segwidth)))
  46.             strtpnt = strtpnt + 2 * segwidth
  47.           END DO
  48.         END IF
  49.       END IF
  50.  
  51.       area = area1 + area2
  52.       ans = area
  53.  
  54.       END !SUBROUTINE  Integrate Function
  55.  
  56.  
  57.  
  58.  
  59.       SUBROUTINE IntegrateVector (y, samper, xl, xh, ans)
  60.       INCLUDE 'STDHDR.FOR'
  61.       REAL y(0:maxv), samper, ans
  62.       INTEGER xl, xh
  63.  
  64.       REAL strtpnt, endpnt, area1, area2, area, segwidth
  65.  
  66.  
  67.       numseg = xh - xl
  68.       segwidth = samper
  69.       endpnt = xh
  70.       area1 = 0.0
  71.       area2 = 0.0
  72.       area = 0.0
  73.       IF (MOD(numseg, 2) .NE. 0) THEN
  74.       !! SIMPSONS 3/8 rule
  75.         area1 = 3.0 / 8.0 * segwidth * (y(endpnt - 3) +
  76.      +          3 * y(endpnt - 2) + 3 * y(endpnt - 1) + y(endpnt))
  77.         endpnt = endpnt - 3
  78.       ELSE
  79.         area1 = 0.0
  80.       END IF
  81.  
  82.       IF (numseg .NE. 3) THEN
  83.       !!SIMPSONS 1/3 rule
  84.         strtpnt = xl
  85.         IF (strtpnt .GE. endpnt) THEN
  86.           area2 = area2 + 1.0 / 3.0 * segwidth *
  87.      +              (y(strtpnt) + 4 * y(strtpnt + 1) + y(strtpnt + 2))
  88.             strtpnt = strtpnt + 2
  89.         ELSE
  90.           DO WHILE (strtpnt .LT. endpnt)
  91.              area2 = area2 + 1.0 / 3.0 * segwidth *
  92.      +              (y(strtpnt) + 4 * y(strtpnt + 1) + y(strtpnt + 2))
  93.             strtpnt = strtpnt + 2
  94.           END DO
  95.         END IF
  96.       END IF
  97.  
  98.       area = area1 + area2
  99.       ans = area
  100.       END !SUBROUTINE
  101.  
  102.       SUBROUTINE RungeKutta(XKnown, YKnown, x, et, h,
  103.      +     HMin, YEst, status)
  104.       REAL XKnown, YKnown, x, et, h, HMin, YEst
  105.       INTEGER status
  106.  
  107.       PARAMETER ( NearZero = 1E-20)
  108.       REAL  xc2,xc3,xc4,xc5,xc6,dc21,dc31,dc32,dc41,dc42,dc43
  109.       REAL  dc51,dc52,dc53,dc54,dc61,dc62,dc63,dc64,dc65
  110.       REAL  yc1,yc3,yc4,yc5,ec1,ec3,ec4,ec5,ec6,d1,d2,d3,d4,d5,d6
  111.       REAL  errr,XEst,EtDiv32
  112.  
  113.       !Calculate the coefficients for the RKF equations
  114.  
  115.       xc2 = .25
  116.       xc3 = 3.0 / 8.0
  117.       xc4 = 12.0 / 13.0
  118.       xc5 = 1.0
  119.       xc6 = .5
  120.       dc21 = .25
  121.       dc31 = 3.0 / 32.0
  122.       dc32 = 9.0 / 32.0
  123.       dc41 = 1932.0 / 2197.0
  124.       dc42 = 7200.0 / 2197.0
  125.       dc43 = 7296.0 / 2197.0
  126.       dc51 = 439.0 / 216.0
  127.       dc52 = 8.0
  128.       dc53 = 3680.0 / 513.0
  129.       dc54 = 845.0 / 4104.0
  130.       dc61 = 8.0/ 27.0
  131.       dc62 = 2.0
  132.       dc63 = 3544.0 / 2565.0
  133.       dc64 = 1859.0 / 4104.0
  134.       dc65 = 11.0 / 40.0
  135.       yc1 = 25.0 / 216.0
  136.       yc3 = 1408.0 / 2565.0
  137.       yc4 = 2197.0 / 4104.0
  138.       yc5 = 0.2
  139.       ec1 = 1.0 / 360.0
  140.       ec3 = 128.0 / 4275.0
  141.       ec4 = 2197.0 / 75240.0
  142.       ec5 = 0.02
  143.       ec6 = 2.0 / 55.0
  144.  
  145.       EtDiv32 = et / 32.0
  146.       XEst = XKnown
  147.       YEst = YKnown
  148.       st = 0
  149.       DO WHILE (st .EQ. 0 )
  150.         d1 = YPrime(XEst, YEst)
  151.         d2 = YPrime(XEst + h * xc2, YEst + h * dc21 * d1)
  152.         d3 = YPrime(XEst + h * xc3, YEst + h * (dc31 * d1 + dc32 * d2))
  153.         d4 = YPrime(XEst + h * xc4, YEst + h *
  154.      +     (dc41 * d1 - dc42 * d2 + dc43 * d3))
  155.         d5 = YPrime(XEst + h * xc5, YEst + h *
  156.      +     (dc51 * d1 - dc52 * d2 + dc53 * d3 - dc54 * d4))
  157.         d6 = YPrime(XEst + h * xc6, YEst + h *
  158.      +     (-dc61 * d1 + dc62 * d2 - dc63 * d3 + dc64 * d4 - dc65*d5))
  159.         errr = h * (ec1 * d1 - ec3 * d3 - ec4 *
  160.      +       d4 + ec5 * d5 + ec6 * d6)
  161.  
  162.         IF (errr .LT. (h * et)) THEN
  163.           XEst = XEst + h
  164.           YEst = YEst + h * (yc1*d1+yc3 * d3 + yc4 * d4 - yc5 * d5)
  165.           IF (ABS(x - XEst) .GT. NearZero) THEN
  166.             IF (errr .LE. (h * EtDiv32))  h = h * 2
  167.             IF (h .GT. (x - XEst))   h = x - XEst
  168.           ELSE
  169.             st = 1
  170.           END IF
  171.         ELSE
  172.           h = h / 2.0
  173.           IF (h .LT. HMin)  st = 2
  174.         END IF
  175.         IF (st .EQ. 2) THEN
  176.           status = 0
  177.         ELSE
  178.           status = 1
  179.         END IF
  180.       END DO !while loop
  181.  
  182.       END !SUBROUTINE {procedure}
  183.  
  184.       FUNCTION YPrime (x, y)
  185.       REAL x, y
  186.         YPrime = EXP(2.0 * x) + y
  187.       END !FUNCTION
  188.  
  189.