- #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