C C PAGE 187-190: NUMERICAL MATHEMATICS AND COMPUTING, CHENEY/KINCAID, 1985 C C FILE: SIMP.FOR C C ADAPTIVE SCHEME FOR SIMPSON'S RULE (SIMP,ASMP,PUSH,POP,FCN) C PARAMETER (IST=5,IFS=16) REAL STACK(IST,3),FSTACK(IFS,3) EXTERNAL FCN DATA EPSI/5.0E-5/, LVMAX/4/ A = 0.0 B = 8.0*ATAN(1.0) CALL ASMP(FCN,A,B,EPSI,LVMAX,STACK,IST,FSTACK,IFS,SUM,IFLAG) PRINT 3,SUM IF(IFLAG .EQ. 0) THEN PRINT 4 ELSE PRINT 5 DO 2 I=1,IFLAG PRINT 6,FSTACK(I,1),FSTACK(I,2),INT(FSTACK(I,3)) 2 CONTINUE END IF 3 FORMAT(//5X,'APPROXIMATE INTEGRAL =',E22.14/) 4 FORMAT(10X,'WITH NO BAD SUBINTERVALS') 5 FORMAT(10X,'WITH BAD SUBINTERVALS:') 6 FORMAT(10X,'[',F10.5,',',F10.5,']',2X,'LEVEL =',I5) STOP END FUNCTION FCN(X) FCN = COS(2.0*X)/EXP(X) RETURN END SUBROUTINE PUSH(A,B,LEVEL,STACK,IST,I) REAL STACK(IST,3) IF(I .LT. IST) THEN I = I+1 STACK(I,1) = A STACK(I,2) = B STACK(I,3) = REAL(LEVEL) RETURN ELSE STOP 'STACK OVERFLOW IN PUSH' END IF END SUBROUTINE POP(A,B,LEVEL,STACK,IST,I) REAL STACK(IST,3) IF(I .GT. 0) THEN A = STACK(I,1) B = STACK(I,2) LEVEL = INT(STACK(I,3)) I = I-1 RETURN ELSE STOP 'STACK UNDERFLOW IN POP' END IF END SUBROUTINE ASMP(FCN,A,B,EPSI,LVMAX,STACK,IST,FSTACK,IFS,SUM,IFLAG) REAL STACK(IST,3),FSTACK(IFS,3) COMMON /BLK/ ITOP,LEVEL,LMAX EXTERNAL FCN ITOP = 0 LEVEL = 0 IFLAG = 0 SUM = 0.0 LMAX = 2**LVMAX IF(IST .GE. LVMAX+1 .AND. IFS .GE. LMAX) THEN CALL PUSH(A,B,LEVEL,STACK,IST,ITOP) DO 2 LOOP=1,2*LMAX-1 IF(ITOP .EQ. 0) RETURN CALL POP(A,B,LEVEL,STACK,IST,ITOP) CALL SIMP(FCN,A,B,EPSI,LVMAX,STACK,IST,FSTACK,IFS,SUM,IFLAG) 2 CONTINUE ELSE PRINT 3 END IF RETURN 3 FORMAT(//5X,'NOT ENOUGH WORKSPACE IN STACK OR FSTACK') END SUBROUTINE SIMP(FCN,A,B,EPSI,LVMAX,STACK,IST,FSTACK,IFS,SUM,IFLAG) REAL STACK(IST,3),FSTACK(IFS,3),X(0:4),F(0:4) COMMON /BLK/ ITOP,LEVEL,LMAX H = B - A H4 = H/4.0 DO 2 I = 0,4 X(I) = A + REAL(I)*H4 F(I) = FCN(X(I)) 2 CONTINUE SUM2 = H*(F(0) + 4.0*F(2) + F(4))/6.0 SUM4 = H*(F(0) + 4.0*F(1) + 2.0*F(2) + 4.0*F(3) + F(4))/12.0 IF(ABS(SUM4 - SUM2) .LE. 15.0*EPSI/REAL(2**LEVEL)) THEN SUM = SUM + SUM4 ELSE IF(LEVEL .LT. LVMAX) THEN LEVEL = LEVEL + 1 CALL PUSH(X(2),X(4),LEVEL,STACK,IST,ITOP) CALL PUSH(X(0),X(2),LEVEL,STACK,IST,ITOP) ELSE CALL PUSH(X(0),X(4),LEVEL,FSTACK,IFS,IFLAG) END IF END IF RETURN END