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 using ! a subroutine. The routines for generating the table ! (getxcoord and getycoord) still need to be subroutines ! because they output a vector. ! ! The trap rule can be programmed as a function because ! it only outputs one double value. ! INTEGER,PARAMETER :: NMAX = 5001 INTEGER :: n REAL(kind=dp) :: a,b,h,total,error,exact REAL(kind=dp) :: xvector(NMAX),yvector(NMAX) ! Declare the trapfun name here. REAL(kind=dp) :: trapfun ! Define problem limits a = -1 b = 4 ! Get the value of n WRITE(*,*) 'Input n' READ(*,*) n ! Call routine to generate x coordinates. This still needs to be ! a subroutine because it outputs a vector CALL getxcoord(n,a,b,xvector) ! Call routine to generate y coordinates ! This also needs to be a subroutine CALL getycoord(n,xvector,yvector) ! Get h from the xvector (can also use h = (b-a)/n here) h = xvector(2)-xvector(1) ! Use an assignment statement to invoke the function. total = trapfun(n,h,yvector) ! 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! FUNCTION trapfun(n,h,y) ! ! 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) :: trapfun INTEGER :: i trapfun = 0 DO i = 1,n+1 IF((i == 1).OR.(i == n+1)) THEN trapfun = trapfun + y(i) ELSE trapfun = trapfun + 2*y(i) ENDIF ENDDO trapfun = trapfun*h/2 RETURN END FUNCTION trapfun