http://www.odu.edu/~agodunov/computing/programs/index.html
Program main
!=============================================
! Integration of a function using Simpson rule
!=============================================
implicit none
double precision f, a, b, integral
integer n, i
double precision, parameter:: pi = 3.1415926
external f
a = 0.0
b = pi
n = 2
write(*,100)
do i=1,16
call simpson(f,a,b,integral,n)
write (*,101) n, integral
n = n*2
end do
100 format(' nint Simpson')
101 format(i9,1pe15.6)
end
Function f(x)
!----------------------------------------
! Function for integration
!----------------------------------------
implicit none
double precision f, x
f = sin(x)
! f = x*cos(10.0*x**2)/(x**2 + 1.0)
return
end
Subroutine simpson(f,a,b,integral,n)
!==========================================================
! Integration of f(x) on [a,b]
! Method: Simpson rule for n intervals
! written by: Alex Godunov (October 2009)
!----------------------------------------------------------
! IN:
! f - Function to integrate (supplied by a user)
! a - Lower limit of integration
! b - Upper limit of integration
! n - number of intervals
! OUT:
! integral - Result of integration
!==========================================================
implicit none
double precision f, a, b, integral,s
double precision h, x
integer nint
integer n, i
! if n is odd we add +1 to make it even
if((n/2)*2.ne.n) n=n+1
! loop over n (number of intervals)
s = 0.0
h = (b-a)/dfloat(n)
do i=2, n-2, 2
x = a+dfloat(i)*h
s = s + 2.0*f(x) + 4.0*f(x+h)
end do
integral = (s + f(a) + f(b) + 4.0*f(a+h))*h/3.0
return
end subroutine simpson
Ada komentar, kritik, saran, atau request?