*DECK DAVINT SUBROUTINE DAVINT (X, Y, N, XLO, XUP, ANS, IERR) C***BEGIN PROLOGUE DAVINT C***PURPOSE Integrate a function tabulated at arbitrarily spaced C abscissas using overlapping parabolas. C***LIBRARY SLATEC C***CATEGORY H2A1B2 C***TYPE DOUBLE PRECISION (AVINT-S, DAVINT-D) C***KEYWORDS INTEGRATION, QUADRATURE, TABULATED DATA C***AUTHOR Jones, R. E., (SNLA) C***DESCRIPTION C C Abstract C DAVINT integrates a function tabulated at arbitrarily spaced C abscissas. The limits of integration need not coincide C with the tabulated abscissas. C C A method of overlapping parabolas fitted to the data is used C provided that there are at least 3 abscissas between the C limits of integration. DAVINT also handles two special cases. C If the limits of integration are equal, DAVINT returns a C result of zero regardless of the number of tabulated values. C If there are only two function values, DAVINT uses the C trapezoid rule. C C Description of Parameters C The user must dimension all arrays appearing in the call list C X(N), Y(N) C C Input-- C X - DOUBLE PRECISION array of abscissas, which must be in C increasing order. C Y - DOUBLE PRECISION array of function values. i.e., C Y(I)=FUNC(X(I)) C N - The integer number of function values supplied. C N .GE. 2 unless XLO = XUP. C XLO - DOUBLE PRECISION lower limit of integration C XUP - DOUBLE PRECISION upper limit of integration. Must have C XLO.LE.XUP C C Output-- C ANS - Double Precision computed approximate value of integral C IERR - A status code C --Normal Code C =1 Means the requested integration was performed. C --Abnormal Codes C =2 Means XUP was less than XLO. C =3 Means the number of X(I) between XLO and XUP C (inclusive) was less than 3 and neither of the two C special cases described in the abstract occurred. C No integration was performed. C =4 Means the restriction X(I+1).GT.X(I) was violated. C =5 Means the number N of function values was .lt. 2. C ANS is set to zero if IERR=2,3,4,or 5. C C DAVINT is documented completely in SC-M-69-335 C Original program from *Numerical Integration* by Davis & Rabinowitz C Adaptation and modifications by Rondall E Jones. C C***REFERENCES R. E. Jones, Approximate integrator of functions C tabulated at arbitrarily spaced abscissas, C Report SC-M-69-335, Sandia Laboratories, 1969. C***ROUTINES CALLED XERMSG C***REVISION HISTORY (YYMMDD) C 690901 DATE WRITTEN C 890831 Modified array declarations. (WRB) C 890831 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE DAVINT C INTEGER I, IERR, INLFT, INRT, ISTART, ISTOP, N DOUBLE PRECISION A, ANS, B, C, CA, CB, CC, FL, FR, R3, RP5, 1 SLOPE, SUM, SYL, SYL2, SYL3, SYU, SYU2, SYU3, TERM1, TERM2, 2 TERM3, X, X1, X12, X13, X2, X23, X3, XLO, XUP, Y DIMENSION X(*),Y(*) C BEGIN BLOCK PERMITTING ...EXITS TO 190 C BEGIN BLOCK PERMITTING ...EXITS TO 180 C***FIRST EXECUTABLE STATEMENT DAVINT IERR = 1 ANS = 0.0D0 IF (XLO .GT. XUP) GO TO 160 IF (XLO .EQ. XUP) GO TO 150 IF (N .GE. 2) GO TO 10 IERR = 5 CALL XERMSG ('SLATEC', 'DAVINT', + 'LESS THAN TWO FUNCTION VALUES WERE SUPPLIED.', + 4, 1) C ...............EXIT GO TO 190 10 CONTINUE DO 20 I = 2, N C ............EXIT IF (X(I) .LE. X(I-1)) GO TO 180 C ...EXIT IF (X(I) .GT. XUP) GO TO 30 20 CONTINUE 30 CONTINUE IF (N .GE. 3) GO TO 40 C C SPECIAL N=2 CASE SLOPE = (Y(2) - Y(1))/(X(2) - X(1)) FL = Y(1) + SLOPE*(XLO - X(1)) FR = Y(2) + SLOPE*(XUP - X(2)) ANS = 0.5D0*(FL + FR)*(XUP - XLO) C ...............EXIT GO TO 190 40 CONTINUE IF (X(N-2) .GE. XLO) GO TO 50 IERR = 3 CALL XERMSG ('SLATEC', 'DAVINT', + 'THERE WERE LESS THAN THREE FUNCTION VALUES ' // + 'BETWEEN THE LIMITS OF INTEGRATION.', 4, 1) C ...............EXIT GO TO 190 50 CONTINUE IF (X(3) .LE. XUP) GO TO 60 IERR = 3 CALL XERMSG ('SLATEC', 'DAVINT', + 'THERE WERE LESS THAN THREE FUNCTION VALUES ' // + 'BETWEEN THE LIMITS OF INTEGRATION.', 4, 1) C ...............EXIT GO TO 190 60 CONTINUE I = 1 70 IF (X(I) .GE. XLO) GO TO 80 I = I + 1 GO TO 70 80 CONTINUE INLFT = I I = N 90 IF (X(I) .LE. XUP) GO TO 100 I = I - 1 GO TO 90 100 CONTINUE INRT = I IF ((INRT - INLFT) .GE. 2) GO TO 110 IERR = 3 CALL XERMSG ('SLATEC', 'DAVINT', + 'THERE WERE LESS THAN THREE FUNCTION VALUES ' // + 'BETWEEN THE LIMITS OF INTEGRATION.', 4, 1) C ...............EXIT GO TO 190 110 CONTINUE ISTART = INLFT IF (INLFT .EQ. 1) ISTART = 2 ISTOP = INRT IF (INRT .EQ. N) ISTOP = N - 1 C R3 = 3.0D0 RP5 = 0.5D0 SUM = 0.0D0 SYL = XLO SYL2 = SYL*SYL SYL3 = SYL2*SYL C DO 140 I = ISTART, ISTOP X1 = X(I-1) X2 = X(I) X3 = X(I+1) X12 = X1 - X2 X13 = X1 - X3 X23 = X2 - X3 TERM1 = Y(I-1)/(X12*X13) TERM2 = -Y(I)/(X12*X23) TERM3 = Y(I+1)/(X13*X23) A = TERM1 + TERM2 + TERM3 B = -(X2 + X3)*TERM1 - (X1 + X3)*TERM2 1 - (X1 + X2)*TERM3 C = X2*X3*TERM1 + X1*X3*TERM2 + X1*X2*TERM3 IF (I .GT. ISTART) GO TO 120 CA = A CB = B CC = C GO TO 130 120 CONTINUE CA = 0.5D0*(A + CA) CB = 0.5D0*(B + CB) CC = 0.5D0*(C + CC) 130 CONTINUE SYU = X2 SYU2 = SYU*SYU SYU3 = SYU2*SYU SUM = SUM + CA*(SYU3 - SYL3)/R3 1 + CB*RP5*(SYU2 - SYL2) + CC*(SYU - SYL) CA = A CB = B CC = C SYL = SYU SYL2 = SYU2 SYL3 = SYU3 140 CONTINUE SYU = XUP ANS = SUM + CA*(SYU**3 - SYL3)/R3 1 + CB*RP5*(SYU**2 - SYL2) + CC*(SYU - SYL) 150 CONTINUE GO TO 170 160 CONTINUE IERR = 2 CALL XERMSG ('SLATEC', 'DAVINT', + 'THE UPPER LIMIT OF INTEGRATION WAS NOT GREATER ' // + 'THAN THE LOWER LIMIT.', 4, 1) 170 CONTINUE C ......EXIT GO TO 190 180 CONTINUE IERR = 4 CALL XERMSG ('SLATEC', 'DAVINT', + 'THE ABSCISSAS WERE NOT STRICTLY INCREASING. MUST HAVE ' // + 'X(I-1) .LT. X(I) FOR ALL I.', 4, 1) 190 CONTINUE RETURN END