SIMPSON

 !Program name: SIMPSON

!Integration of a function from Simpson’s rule

!Integral of the function exp(x) from 0 to 1 interval

!the exact integration is exp(1)-1=1.1718282

!h=(b-a)/2

!x1=a+h

! S=h/3[f(b)+4f(x1)+f(a)}

! n: number of interval (even number)

!S=integral

 

Real :: fx, exact, integral, a, b

Integer :: i, n

exact=exp(1.)-1

write (*,*) "Exact integration of exp(x) from 0 to 1= ", exact

a=0

b=1

n=2

call simpson(a,b,h,fx1,n,integral)

print *, n, h, fx1, integral

stop

end

 

subroutine simpson(a,b,h,fx1,n,integral)

real :: fa, fb, fx1, integral, h, a, b

integer :: n, i

fa=exp(a)

fb=exp(b)

h=(b-a)/dfloat(n)

fx1=exp(a+h)

integral = (fa+fb+4.0*fx1)*h/3.0

return

end subroutine simpson