Search code examples
fortrannumericalcalculus

Problems with the Trapezoidal Rule


I'm having some troubles to calcule the integral of e^x inside and interval [b.a] using fortran.

I think I'm doing something wrong in the funcion calls. Thanks for helping me.

program trapezium
  implicit none

    integer :: i, n, b, a
    real :: sumation, mean, deltax, f(i), integral


 ! The value of the integral using the trapezium rule can be found using
 ! integral = (b - a)*((f(a) +f(b))/2 + sumation_1_n-1 )/n 

write(*,*) "type the limits b, a and the number of intervals"
     read *, b, a, n

    deltax = (b - a)/n
        mean = (f(a) + f(b))/2
sumation = 0

do i = 1, n-1  
    sumation = sumation + f(i)
end do


      integral = deltax*(mean + sumation) 
  write (*,*) "the value of the integral using the trapezoidal method is", integral

     end program 

function f(x)
  real :: f(x) 
  integer :: x

      f(x) = EXP(x)

end function

Solution

  • There are a couple of issues with your code:

    • f is a function, but at the same time you define an array f(i)
    • When defining an array of fixed size, the size has to be known at compile time. So real :: f(i) is only valid for a constant i
    • exp() expects a real variable, not an integer
    • Integer arithmetic might lead to unexpected results: 1/2 = 0 and not 0.5!

    What about (This does not try to fix the maths, though - see my comment):

    module functions
    contains
      function f(x)
        implicit none
        real :: f
        integer,intent(in) :: x
    
        f = EXP(real(x))
    
      end function
    end module
    
    program trapezium
      use functions
      implicit none
    
      integer :: i, n, b, a
      real :: sumation, mean, deltax, integral
    
    
      ! The value of the integral using the trapezium rule can be found using
      ! integral = (b - a)*((f(a) +f(b))/2 + sumation_1_n-1 )/n 
    
      write(*,*) "type the limits b, a and the number of intervals"
      read *, b, a, n
    
      deltax = real(b - a)/real(n)
      mean = (f(a) + f(b))/2
      sumation = 0
    
      do i = 1, n-1  
        sumation = sumation + f(i)
      end do
    
    
      integral = deltax*(mean + sumation) 
      write (*,*) "the value of the integral using the trapezoidal method is", integral
    
    end program 
    

    Note that the use of modules enables the compiler to check for the arguments of the function. Additionally, you do not need to define the return value of the function in the main program.