a fortran environment

[danbaron]

Here is a shorter version of the previous program.

I don't think Simpson's Rule is the fastest way to numerically integrate, but, you can get pretty good accuracy if you divide the interval into tiny sub-intervals.

Below the interval is divided into ten million equal sub-intervals, and the result is accurate to 1 part in 3,141,592,653,589 - approximately one part in three trillion.

' code ----------------------------------------------------------------------------------------------------------------

module integrate

contains

double precision function xsquaredcosxdiv4(x)
double precision::x
xsquaredcosxdiv4=x*x*cos(x)/4
end function

double precision function simpsonintegrate(x1,x2,ndivisions,f)
double precision::x1,x2,f
integer::ndivisions
double precision::division,sum
integer::i
division=(x2-x1)/ndivisions
sum=(f(x1)+f(x2))/2
do i= 1,ndivisions-1
sum=sum+f(x1+division*i)
enddo
simpsonintegrate=division*sum
end function

end module

!-----------------------------------------------------------------------------------------------------------

program test
use integrate
implicit none
double precision::pi
pi=4*atan(1d0)
print*,"pi               = ",pi
print*,"integral approx. = ",simpsonintegrate(0d0,2*pi,10000000,xsquaredcosxdiv4)
end program

' output --------------------------------------------------------------------------------------------------------------

 pi               =    3.1415926535897931
 integral approx. =    3.1415926535892416