next up previous
Next: Error estimate Up: Integration Previous: Integration

integrate.f

c      PROGRAM TO INTEGRATE  f(x) OVER [a,b]
        program integrate
       implicit none
       real*8 a,b
       integer nsubs,nxs
       integer i
       real*8 x,dx
       real*8 area
       dimension x(101)

c      INTERVAL [a,b] IS DIVIDED BY nxs EQUALLY-SPACED X-VALUES
c      INTO nsubs SUBINTERVALS, EACH OF LENGTH dx
c      x IS A VECTOR CONTAINING THESE X-VALUES
c      THE SUBROUTINE find_int IS CALLED TO DO THE AREA COMPUTATION
c      AND THE RESULT IS PRINTED USING print

       open(unit=1,file='int.dat')
        read(1,2000)a
        read(1,2000)b
        read(1,2001)nsubs
       close(1)

       if (nsubs .gt. 100) then
         write(*,*)"Only allowed 100 subintervals...try again"
         go to 999
       endif

       nxs=nsubs+1

       dx=1.d0/nsubs
       do i=1,nxs
         x(i)=(i-1)*dx
       enddo

       call find_int(x,dx,nxs,area)

       call print(nsubs,area)

2000   format('',f12.7)
2001   format('',I3)
999    stop
       end


c     7**************************************************************72
       subroutine find_int(x,dx,nxs,area)
       implicit none
       integer nxs
       real*8 x,dx
       real*8 ff
       real*8 area
       integer i
       dimension x(1)
c      THIS SUBROUTINE DOES THE WORK IN FINDING THE AREA
c      UNDER f OVER [a,b].  THE INTEGRAL IS COMPUTED USING
c      LEFT ENDPOINTS.
c      x IS A VECTOR CONTAINING nxs X-VALUES

       area=0.d0
       do i=1,nxs-1
         area=area+dx*ff(x(i))
       enddo

       return
       end

c     7**************************************************************72
       subroutine print(nsubs,area)
       implicit none
       integer nsubs
       real*8 area
c      THIS SUBROUTINE WRITES THE OUTPUT TO A FILE

       open(unit=1,file='integral_out')
       write(1,2000)nsubs,area

2000   format('','COMPUTING THE INTEGRAL WITH ',I4,
     +  ' POINTS, THE AREA IS ',f12.7)

       return
       end

c     7**************************************************************72
       real*8 function ff(x)
c      THIS FUNCTION SUBROUTINE COMPUTES THE VALUE OF f AT A POINT x
       implicit none
       real*8 x

       ff = x*x

       return
       end

This routine reads the file int.dat. This file needs to look like

0.00000000
1.00000000
100

giving the interval endpoints and the number of subintervals.

To get better accuracy, one can use the trapezodial rule. Instead of using just the left or right endpoint value, average these values. Geometrically, this means that instead of using rectangles, use the trapezoid whose top is the line between the values of f at the left and right. As an assignment, change the code to a trapezodial algorithm.



Bruce Pitman
Wed Oct 11 12:23:54 EDT 1995