! compile 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