How to Find the Definite Integral of a Function Using the Trapezoidal Rule

The trapezoidal rule is one of the easiest ways to find a definite integral. It belongs to the class of formulas called Newton-Cotes formulas which evaluate the integrand at equally spaced points. The trapezoidal rule is extremely accurate for periodic functions when they are integrated over their periods. Let us see how to evaluate a definite integral using the trapezoidal rule.

• We read the limits of integration, a and b at the time of execution.
• Suppose we want the integral for many values of n (the number of divisions). We can call a subroutine to integrate the function for each value of n using a DO loop. We calculate h=(b-a)/n, the size of the interval, for each value of n. In general, the higher the value of n, closer the upper line of the trapezium to the actual function, and the better the results. (See the image)
• We need to define the function, f(x) we wish to integrate. We will need this in the subroutine.
•  We create a subroutine to give the value of the integral.
• We calculate sum1=[f(a)+f(b)]/2. Since sum1 is of a fixed value, it can also be calculated outside the subroutine. But, let's not get there.
• For k=1 to (n-1), we calculate temp=f(a+k*h) and sum up the values in the variable temp itself.
• We add sum1 and temp and multiply the result by the length of the interval, h using the expression, int_trap=h*(sum1+temp) where int_trap is the value of the integral.
• Optional: To check the results (assuming the function used is integrable), we define another function, g(x) which is the integral of f(x). The value of int_trap and another variable check=g(b)-g(a) should be the same. We check the absolute value of (int_trap-check). If it is <<0 we display that the results match. Else we ask to make adjustments or increase the size of n for higher accuracy.

Let us look at the derivation of the formula for trapezoidal rule used in the program.

Algorithm Used in the Program for Integrating sin(x) from 0 to 1. We have used the values of angles in radians. So 1 means 1 Radian. Note that the answer is the same as -[cos(1)-cos(0)]. This is because the integral of sin(x) is cos(x).

Now for the Fortran Program,

!To integrate a function over an interval using Trapezoidal ruleprogram trapezoidal
    implicit none
    integer n
    real a,b,h,int_trap,check,g,diff !The function 'g' has to be declared here
   
    write(*,*) "Enter the values of a and b to specify the limits of integration"
    write(*,*) "Enter values in radians" !Optional. Only for trigonometric functions
    read (*,*) a,b
    write(*,*) "a= ",a,"b= ",b
   
    write(*,*) "The given function is f(x)=sin(x)" !Change this if you change the function
    write(*,*)
    write(*,10) "n", "int-trap" !Headings for the result
    10 format (3x,a10,3x,a10)
    write(*,*)
   
    do n=5,200,5
    h=(b-a)/n
    call trap(a,b,n,h,int_trap) !To call the integral value
    write (*,*) n,int_trap !To write the integral for each 'n'
    end do
   
    check=g(b)-g(a) !Optional. Skip to the subroutine if you do not wish to verify results
    write(*,*)
    write(*,*) "The value using the integral function, [-cos(x)] is= ",check !Change the function in the character string if you change f(x)
   
    n=200
     call trap(a,b,n,h,int_trap)
    diff=(abs(check-int_trap))
    if (diff.lt.0.00001) then
          write (*,*) "For n= ",n," The difference in the two answers is =",diff,". The answer has been verified."
          else
        write(*,*) "For n= ",n," The difference in the two answers is =",diff,". Trapezoidal rule fails for this integral."
        write(*,*) "Try increasing the size of n for more accurate results." !To prompt user to increase 'n' for better results
    end if
   
end program trapezoidal

   
    !To calculate the integral
    subroutine trap(a,b,n,h,int_trap)
        implicit none
           real a,b,h,sum1,int_trap,temp,f !The function f has to be declared here
        integer n,k
   
        sum1=0.0
        temp=0.0
        sum1=(f(a)+f(b))/2
   
        do k=1,n-1
            temp=temp+f(a+(k*h))
        end do
   
        int_trap=h*(sum1+temp)
    end subroutine trap
   

    !To define the function to integrate
    real function f(x)
        implicit none
          real::x
      
          f=sin(x) !Change this to the function of your choice
          return
    end function f
       

    !To define the integral of the above function. Optional:Use only if you want to verify results
    real function g(x)
        implicit none
          real::x
          g=-cos(x) !Change this depending on f(x)
          return
    end function g

Sample Output for the Trapezoidal Rule. We have used the function sin(x) in the range 0 to 1 radians.