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.