! Pi calculator, by Reginald Ford ! James Madison University ! PLEDGE: The JMU Honor Code was obeyed during this work !---------------------------------------------------------------------- ! The following is a module for simpson's rule along with ! a function, f(x), which is the semicircle function. ! Simpson's rule is used to integrate over the domain of the function, ! and this provides pi*r^2/2 . After multiplying by (2/r^2), the approximation of ! pi is printed. Several approximations are done to see the convergence. !---------------------------------------------------------------------- module simpsons_rule ! r2 should be a perfect square so that there is ! no error in computing r, which is used in the bounds of integration ! This value seems to have no effect on the accuracy of the calculations real (kind=8)::r2=2**30 contains ! The semicircle function for radius sqrt(r2) real (kind = 8) function f(x) implicit none real (kind=8),intent(in)::x f=sqrt(r2-x**2.) ! sqrt(x) provides nan for x<0 ! So this is just in case if(isnan(f)) then f=0.0 end if return end function f ! This is the simple simpson's rule, which is basically ! quadratic interpolation and integration from a to b ! using (a+b)/2 as the midpoint real (kind=8) function simple(a,b) implicit none real (kind=8),intent(in)::a,b simple=((b-a)/6.)*(f(a)+4*f((a+b)/2.)+f(b)) return end function simple ! This will integrate f(x) from a to b using simpson's rule ! from a to b with 'accuracy' controlling the number of nodes used real (kind=8) function composite(a,b,accuracy) implicit none real (kind=8) , intent(in)::a,b integer, intent(in)::accuracy integer::n,i real (kind=8)::h,integral integral=0.0 n=accuracy ! the number of nodes needs to be odd if(mod(n,2)==0) then n=n+1 end if h=(b-a)/n do i=0,n-1 integral=integral+simple(a+i*h,a+(i+1)*h) end do composite=integral return end function composite end module simpsons_rule ! MAIN PROGRAM STARTS HERE ! program pi_calc use simpsons_rule implicit none integer::max_t=18 integer::t,n real (kind=8)::approx,r r=sqrt(r2) write(*,*)"Approximations of pi:" do t=1,max_t,2 ! This allows us to see how it converges, with a large range of n covered n=3**t approx=composite((-1.)*r,r,n)*(2./r2) write(*,*)"n=",n,":",approx," error:",acos(-1.)-approx end do ! The program takes too long to converge to epsilon(double) precision on my machine write(*,*)"epsilon=",epsilon(r) end program pi_calc