module simpson double precision eps end module simpson Program main !------------------------------------------------------------ ! Integration of a function f(x,y) using Simpson rule twice ! with doubling number of intervals and recursive calling !------------------------------------------------------------ use simpson implicit none double precision a, b, f, simpson2D, integral external f a = 0.0 b = 1.0 eps = 1.0e-7 write(*,100) integral = simpson2D(f,a,b,eps) write(*,103) integral 100 format(' 2D integration + Simpson rule') 103 format(' integral = ',f12.8) stop End Function f(x) !---------------------------------------- ! Function for integration over x !---------------------------------------- use simpson implicit none double precision f, g, x, x1 double precision c, d, simpson2D common/fxy/x1 external g c = 0.0 d = sin(x) x1 = x ! x1 passes the value of x to g(y) f = simpson2D(g,c,d,eps) return end Function g(y) !---------------------------------------- ! Function for integration over y !---------------------------------------- implicit none double precision g, y, x common/fxy/x ! g = (x + y)**2 g = x*x/(y*y+2.0) return end recursive function simpson2D(f,a,b,eps) !========================================================== ! 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 !========================================================== implicit none double precision f, a, b, eps, simpson2D double precision sn, s2n, h, x 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 simpson2D = s2n exit end if sn = s2n n = n*2 end do return end function simpson2D