PROGRAM D4R3 C Driver for routine QSIMP PARAMETER(PIO2=1.5707963) A=0.0 B=PIO2 WRITE(*,'(1X,A)') 'Integral of FUNC computed with QSIMP' WRITE(*,'(1X,A,F10.6)') 'Actual value of integral is', * FINT(B)-FINT(A) CALL QSIMP(A,B,S) WRITE(*,'(1X,A,F10.6)') 'Result from routine QSIMP is',S END FUNCTION FUNC(X) FUNC=(X**2)*(X**2-2.0)*SIN(X) END FUNCTION FINT(X) C Integral of FUNC FINT=4.0*X*((X**2)-7.0)*SIN(X) * -((X**4)-14.0*(X**2)+28.0)*COS(X) END C SUBROUTINE QSIMP(A,B,S) PARAMETER (EPS=0.5E-5, JMAX=200) OST=-1.E30 OS= -1.E30 DO 11 J=1,JMAX CALL TRAPZD(A,B,ST,J) S=(4.*ST-OST)/3. IF (ABS(S-OS).LT.EPS) RETURN OS=S OST=ST 11 CONTINUE PAUSE 'Too many steps.' END C SUBROUTINE TRAPZD(A,B,S,N) IF (N.EQ.1) THEN S=0.5*(B-A)*(FUNC(a)+FUNC(b)) ELSE DEL=(B-A)/n SUM=0.5*(func(a)+func(b)) x=a+del DO 11 J=2,n SUM=SUM+FUNC(x) X=X+DEL 11 CONTINUE S=del*sum ENDIF RETURN end