1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Homework Help: IFORT Coding Problem (Simpson's Rule)

  1. May 11, 2009 #1
    1. The problem statement, all variables and given/known data

    For a function f(x) over the interval [a,b], simpson's rule approximates the defenite integral:

    [tex]\int[/tex]f(x) dx = h/3*[f0 + 4(f1 + f3+ ...+ f2n-1) + 2(f2 + f4 + ... + f2n-2) + f2n]

    where h = (b-a)/n
    fi is f(a + i*h)

    Given an error epsilon, I need to write an executable program using IFORT that will find this integral using the simpson's method for f(x) = exp(-x)*sin(x) on [0,10] which converges if
    new integral - old integral < new integral*eps
    obvioulsy a loop will be in this program. What you have to do is cut the h in half every time to get a better approximation. And each time you have to change the odd fi to the even fi, and vice versa.

    2. Relevant equations

    see the above equation

    3. The attempt at a solution

    Here is my program.

    PROGRAM Simpson_Rule
    REAL(8) :: eps, a, b, integral, f
    eps = 1.0d-10
    a = 0.0d0
    b = 10.0d0
    CALL Simpson_Integrate(f, a, b, eps, integral)
    PRINT *, integral
    END PROGRAM Simpson_Rule

    SUBROUTINE Simpson_Integrate(f, a, b, eps, integral)
    REAL(8) :: a, b, integral, f, eps, h, x, sum, sumeven, sumodd,
    & integral_old, end_points
    INTEGER :: i, n_sum_add, ne
    h = (b-a)/2.0d0
    sum = f(a) + f(b)
    end_points = sum/3.0d0
    x = (a+b)/2.0d0
    sumodd = f(x)
    integral_old = h*end_points + 4.0d0*sumodd*h/3.0d0
    n_sum_add = 2
    do ne = 1, 40
    sumeven = sumeven + sumodd
    sumodd = 0.0d0
    x = a + h/2.0d0
    do i = 1, n_sum_add
    sumodd = f(x)
    x = x+h
    end do
    h = h/2.0d0
    integral = h*end_points + (4.0d0*sumodd + 2.0d0*sumeven)*h/3.0d0
    print*, 2**ne, integral
    if (abs(integral-integral_old) < abs(integral)*eps) return
    integral_old = integral
    n_sum_add = n_sum_add*2
    end do
    print*, "Integration has not converged"
    END SUBROUTINE Simpson_Integrate

    real(8) function f(x)
    REAL(8) :: x
    f = exp(-x)*sin(x)
    end function f

    When I execute, i don't get sensible results. Any help is appreciated
  2. jcsd
  3. May 11, 2009 #2


    User Avatar
    Science Advisor

    The first thing I noticed is that fortran might not like having a function named the same as a variable.

    I was confused when you wrote f(a) + f(b); I originally thought that f was an array. Try changing the name of that function.
  4. May 11, 2009 #3
    I don't believe that is the problem. I've done this same type problem with Trapezoid rule and it works just fine with the same function defined as I have done in the program shown above.
  5. May 11, 2009 #4
    Ok I actually solved this problem on my own. In the loop I need:

    sumodd = f(x) + sumodd

    and then all is swell.

Share this great discussion with others via Reddit, Google+, Twitter, or Facebook