CODEWITHSUNDEEP
fortrannumericalcalculus
user2803219
user2803219

Reputation: 55

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

Upvotes: 2

Views: 8453

Answers (1)

Alexander Vogt
Alexander Vogt

Reputation: 18098

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.

Upvotes: 2

Related Questions