!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Composite Trapezoid Integration ! ! with Romberg Extrapolation ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Consider a function f(x) bounded by [a,b]. !! !! The function is broken up into subintervals !! !! and each subinterval is evaluated with the !! !! trapezoid rule. A Richardson extrapolation !! !! is then used to model the behavior of the !! !! partial sums. Richardson speeds up the !! !! convergence of the integration and greatly !! !! increases the efficiency of the trapezoid !! !! rule. !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! compiled with gfortran -fdefault-real-8 Composite_Trapezoid_Romberg.f95 !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ABW 2012 - 2015 ! !!!!!!!!!!!!!!!!!!!!! program comptrap implicit none integer :: n,k,i,j real(kind=8),parameter :: answ = 0.0887552844352566d0 ! x^2*sin(x) from 0 to pi/4 ! real(kind=8),parameter :: answ = 0.1922590000000000d0 ! x^2*ln(x) from 1 to 1.5 ! real(kind=8),parameter :: answ = 2.5886286325071759d0 ! exp(3*x)*sin(2*x) from 0 to pi/4 ! real(kind=8),parameter :: answ = -0.7339690000000000d0 ! 2*x/(x^2 - 4) from 1 to 1.6 real(kind=8),parameter :: pi = 3.1415926535897932d0 ! Exact value of pi to 16 digits real(kind=8),parameter :: a = 0.0d0,b = pi/4.0d0 ! Bounds of integration [a,b] real(kind=8),dimension(0:3,0:3) :: R,sigR real(kind=8) :: hk,sumf R(0,0) = ((b - a)/2.0d0)*(f(a) + f(b)) do k = 1, 3 hk = (b-a)/(2.0d0**k) sumf = 0.0d0 do i = 1,2**(k - 1) sumf = sumf + f(a + hk*(2*i - 1)) enddo R(k,0) = 0.5d0*R(k-1,0) + hk*sumf enddo do k = 1, 3 do j = 1,k R(k,j) = (1.0d0/(4.0d0**j - 1.0d0))*((4.0d0**j)*(R(k,j-1)) - R(k-1,j-1)) enddo enddo ! Calculate error in Romberg extrapolation do k = 0,3 do j = 0,3 sigR(k,j) = abs(R(k,j) - answ) enddo enddo write(*,101) 'Function to integrate: f(x) = x^2*sin(x) over [',a,', ',b,']' 101 format (A48,F3.1,A2,F7.5,A1) write(*,*) write(*,"(A14)") 'Approximation:' write(*,*) ' 0',' 1',' 2',' 3' write(*,"(A,F12.9,$)") '0 ',R(0,0) write(*,*) write(*,"(A,F12.9,F12.9)") '1 ',R(1,0),R(1,1) write(*,"(A,F12.9,F12.9,F12.9)") '2 ',R(2,0),R(2,1),R(2,2) write(*,"(A,F12.9,F12.9,F12.9,F12.9)") '3 ',R(3,0),R(3,1),R(3,2),R(3,3) write(*,*) write(*,"(A,F12.9)") 'Expected Value:',answ write(*,*) write(*,"(A19,F20.16)") '|R(3,3) - Itrue| =',abs(R(3,3) - answ) write(*,*) write(*,"(A21)") 'Error in Integration:' write(*,*) ' 0',' 1',' 2',' 3' write(*,"(A,F12.9,$)") '0 ',sigR(0,0) write(*,*) write(*,"(A,F12.9,F12.9)") '1 ',sigR(1,0),sigR(1,1) write(*,"(A,F12.9,F12.9,F12.9)") '2 ',sigR(2,0),sigR(2,1),sigR(2,2) write(*,"(A,F12.9,F12.9,F12.9,F12.9)") '3 ',sigR(3,0),sigR(3,1),sigR(3,2),sigR(3,3) contains ! Function to integrate function f(x) implicit none real(kind=8) :: f, x f = x*x*sin(x) ! f = x*x*log(x) ! f = exp(x*3.0d0)*sin(x*2.0d0) ! f = (x*2.0d0)/(x*x - 4.0d0) return end function f end program comptrap