Program main !------------------------------------------------------------------ ! Integration of a function using Simpson rule ! with doubling number of intervals !------------------------------------------------------------------ implicit none double precision a, b, f, integral, eps external f integer nint double precision, parameter :: pi = 3.1415926 a = 0.0 b = pi eps = 1.0e-8 call simpson2(f,a,b,eps,integral,nint) write(*,103) nint, integral 103 format(' intervals = ',i8,/,' integral = ',f15.9) stop 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 simpson2(f,a,b,eps,integral,nint) !========================================================== ! Integration of f(x) on [a,b] ! Method: Simpson rule with doubling number of intervals ! till error = coeff*|I_n - I_2n| < eps ! 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 ! eps - tolerance ! OUT: ! integral - Result of integration ! nint - number of intervals to achieve accuracy !========================================================== implicit none double precision f, a, b, eps, integral double precision sn, s2n, h, x integer nint double precision, parameter :: coeff = 1.0/15.0 ! error estimate coeff integer, parameter :: nmax=1048576 ! max number of intervals integer n, i ! evaluate integral for 2 intervals (three points) h = (b-a)/2.0 sn = (1.0/3.0)*h*(f(a)+4.0*f(a+h)+f(b)) ! loop over number of intervals (starting from 4 intervals) n=4 do while (n <= nmax) s2n = 0.0 h = (b-a)/dfloat(n) do i=2, n-2, 2 x = a+dfloat(i)*h s2n = s2n + 2.0*f(x) + 4.0*f(x+h) end do s2n = (s2n + f(a) + f(b) + 4.0*f(a+h))*h/3.0 if(coeff*abs(s2n-sn) <= eps) then integral = s2n + coeff*(s2n-sn) nint = n exit end if sn = s2n n = n*2 end do return end subroutine simpson2