Program main !------------------------------------------------------------------ ! Integration of a function using adaptive procedure ! with Gauss 8 and 16 points !------------------------------------------------------------------ implicit none double precision gauss_a, a, b, f, integral, eps external f integer nint double precision, parameter :: pi = 3.1415926 a = 0.0 b = 2.0*pi eps = 1.0e-7 integral = gauss_a(f,a,b,eps) write(*,103) integral 103 format(' integral = ',f12.8) 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) f = x*sin(30.0*x)*cos(50.0*x) return end recursive function gauss_a(f,a,b,eps) !========================================================== ! Integration of f(x) on [a,b] ! Method: Gauss quadratures with adaptive integration ! for left/right intervals till error = |I_16 - I_8| < 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 (should not be less than 1.0e-8) ! OUT: ! gauss_a - Result of integration !========================================================== implicit none double precision gauss_a, f, a, b, eps, gauss8, gauss16 double precision s1, s2, h, ax, bx, sum integer i external f if(eps <= 1.0e-8) eps = 1.0e-8 h = (b-a)/2.0 sum = 0.0 do i=1,2 ax = a + h*dfloat(i-1) bx = ax + h s1 = gauss8(f,ax,bx) s2 = gauss16(f,ax,bx) if(abs(s2-s1)<= eps .and. abs(s2-s1)/abs(s2+s1)<= eps) then sum = sum + s2 else sum = sum + gauss_a(f,ax,bx,eps) end if end do gauss_a = sum return end function gauss_a Function gauss8(f,a,b) !========================================================== ! Integration of f(x) on [a,b] ! Method: Gauss 8 points ! 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 ! OUT: ! gauss8 - Result of integration !========================================================== implicit none integer, parameter :: n=4 double precision gauss8, f, a, b double precision ti(n), ci(n) data ti/0.1834346424, 0.5255324099, 0.7966664774, 0.9602898564/ data ci/0.3626837833, 0.3137066458, 0.2223810344, 0.1012285362/ double precision r, m, c integer i r = 0.0; m = (b-a)/2.0; c = (b+a)/2.0; do i = 1,n r = r + ci(i)*(f(m*(-1.0)*ti(i) + c) + f(m*ti(i) + c)) end do gauss8 = r*m return end function gauss8 Function gauss16(f,a,b) !========================================================== ! Integration of f(x) on [a,b] ! Method: Gauss 16 points ! 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 ! OUT: ! gauss16 - Result of integration !========================================================== implicit none integer, parameter :: n=8 double precision gauss16, f, a, b double precision ti(n), ci(n) data ti/0.0950125098, 0.2816035507, 0.4580167776, 0.6178762444, & 0.7554044083, 0.8656312023, 0.9445750230, 0.9894009349/ data ci/0.1894506104, 0.1826034150, 0.1691565193, 0.1495959888, & 0.1246289712, 0.0951585116, 0.0622535239, 0.0271524594/ double precision r, m, c integer i r = 0.0; m = (b-a)/2.0; c = (b+a)/2.0; do i = 1,n r = r + ci(i)*(f(m*(-1.0)*ti(i) + c) + f(m*ti(i) + c)) end do gauss16 = r*m return end function gauss16