From Education

Tom's numerical integration of x^2 in Fortran90

! complile with a command like "gfortran integrate.f90"
! run with the command "./a.out 1000 0 1"
!     to integrate f(x)=x^2 using 1000 steps from 0 to 1
!     the theoretical answer is 1/3

! This is the function "f" being integrated
FUNCTION f(x)
    REAL(8) :: f
    REAL(8) x
    f = x*x
    RETURN
END FUNCTION f

! Function "init" calculates deltaX and initializes yVals
SUBROUTINE init(numSeg, startX, endX, deltaX, yVals)
    INTEGER numSeg, i
    REAL(8) startX, endX, deltaX,  yVals(0:numSeg)
	REAL(8) x, rangeX
    rangeX = endX - startX
    deltaX = rangeX/numSeg
    DO i=0, numSeg
        x = startX + i*rangeX/numSeg
        yVals(i) = f(x)
    END DO
END SUBROUTINE init

!  Add the area of all the trapazoids together, to estimate the integral
FUNCTION    integrate(yVals, numSeg, deltaX)
    REAL(8) :: integrate
    INTEGER numSeg, i
    REAL(8) deltaX, yVals(0:numSeg)
    integrate = 0
    DO i=0, numSeg
        ! WRITE (*, '(F23.18)') area
        integrate = integrate + (yVals(i) + yVals(i+1))*.5*deltaX
    END DO
    RETURN
END FUNCTION integrate

PROGRAM integrateFx
    IMPLICIT NONE
    INTEGER numSeg
	REAL(8) startX, endX, deltaX, integrate
    REAL(8), ALLOCATABLE :: yVals(:)
    CHARACTER *100 BUFFER

!  Get number of segments, starting and ending values for integral
    IF (COMMAND_ARGUMENT_COUNT() < 3) STOP
    CALL GET_COMMAND_ARGUMENT(1, BUFFER)
    READ (BUFFER, *) numSeg
    CALL GET_COMMAND_ARGUMENT(2, BUFFER)
    READ (BUFFER, *) startX
    CALL GET_COMMAND_ARGUMENT(3, BUFFER)
    READ (BUFFER, *) endX
    IF ((numSeg<1) .OR. (startX >= endX)) STOP

!  Allocate space for function values and calculate them
    ALLOCATE(yVals(0:numSeg))
    CALL init(numSeg, startX, endX, deltaX, yVals)

!  Display result
    WRITE (*,'(a28,1x,f6.2,1x,a4,1x,f6.2,1x,a2,1x,g23.18)') &
        'Area under curve y=x^2 from', startX, 'to', endX, 'is', &
        integrate(yVals, numSeg, deltaX)
END PROGRAM