MODULE const INTEGER, PARAMETER :: dp = KIND(1.0d0), sp = KIND(1.0e0) REAL(KIND=dp), PARAMETER :: dpi = 3.141592653589793d0 END MODULE const ! PROGRAM prog USE const IMPLICIT NONE ! ! Main program for implementing the trap rule in ! subroutines. INTEGER,PARAMETER :: NMAX = 5001 INTEGER :: n REAL(kind=dp) :: a,b,h,total,error,exact REAL(kind=dp) :: xvector(NMAX),yvector(NMAX) ! Define problem limits a = -1 b = 4 ! Get the value of n WRITE(*,*) 'Input n' READ(*,*) n ! Call routine to generate x coordinates CALL getxcoord(n,a,b,xvector) ! Call routine to generate y coordinates CALL getycoord(n,xvector,yvector) ! Get h from the xvector (can also use h = (b-a)/n here) h = xvector(2)-xvector(1) ! Do the trap rule CALL traprule(n,h,yvector,total) ! Cmopare approximate value with the exact value exact = -5/EXP(4.0d0) + COS(1.0d0)-COS(4.0d0) error = (total-exact)/exact WRITE(*,*) 'total = ',total WRITE(*,*) 'exact = ',exact WRITE(*,*) 'error = ',error END PROGRAM prog !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE getxcoord(n,a,b,x) ! ! Subroutine to generate a set of equally spaced ! x coordinates given the interval limits [a,b] and ! the number of subdivisions n ! USE const IMPLICIT NONE INTEGER,INTENT(IN) :: n REAL(kind=dp),INTENT(IN) :: a,b REAL(kind=dp),INTENT(OUT) :: x(n+1) INTEGER :: i REAL(kind=dp) :: h h = (b-a)/n DO i = 1,n+1 x(i) = a +(i-1)*h ENDDO RETURN END SUBROUTINE getxcoord !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE getycoord(n,x,y) ! ! Subroutine to get a set of y coordinates that ! go with a set of x coordinates ! USE const IMPLICIT NONE INTEGER,INTENT(IN) :: n REAL(kind=dp),INTENT(IN) :: x(n+1) REAL(kind=dp),INTENT(OUT) :: y(n+1) INTEGER :: i DO i = 1,n+1 y(i) = x(i)*EXP(-x(i)) + SIN(x(i)) ENDDO RETURN END SUBROUTINE getycoord !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE traprule(n,h,y,total) ! ! Subroutine to implement the trapezoidal approximation ! to a definite integral. ! USE const IMPLICIT NONE INTEGER,INTENT(IN) :: n REAL(kind=dp),INTENT(IN) :: h, y(n+1) REAL(kind=dp),INTENT(OUT) :: total INTEGER :: i total = 0 DO i = 1,n+1 IF((i == 1).OR.(i == n+1)) THEN total = total + y(i) ELSE total = total + 2*y(i) ENDIF ENDDO total = total*h/2 RETURN END SUBROUTINE traprule