- #1
Physics_Math
- 13
- 0
Homework Statement
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.
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