IFORT Coding Problem (Simpson's Rule)

[Physics_Math]

Homework Statement

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


$$\int$$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.

Homework Equations

see the above equation

The Attempt at a Solution

Here is my program.

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

SUBROUTINE Simpson_Integrate(f, a, b, eps, integral)
IMPLICIT NONE
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"
return
END SUBROUTINE Simpson_Integrate

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

When I execute, i don't get sensible results. Any help is appreciated

 

[minger]

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.

 

[Physics_Math]

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.

 

[Physics_Math]

Ok I actually solved this problem on my own. In the loop I need:

sumodd = f(x) + sumodd

and then all is swell.

:D