!======================================================================= ! ch07_improper_infinite.f90 ! ! Improper integral with an infinite upper limit. ! ! The program evaluates ! ! infinity ! / ! | exp(-x) sin(x) dx = 1/2 ! / ! 0 ! ! by transforming the semi-infinite interval [0,infinity) to the ! finite interval [0,1): ! ! x = t/(1-t), dx = dt/(1-t)^2 . ! ! The transformed integral is ! ! 1 ! / ! | exp[-t/(1-t)] sin[t/(1-t)] / (1-t)^2 dt . ! / ! 0 ! ! Gauss-Legendre quadrature is then applied on [0,1]. The quadrature ! nodes do not include the endpoints, so the transformed function is ! never evaluated at t = 1. ! ! This example follows Section 7.5.1 of the book: ! "Improper integrals: type 1, infinite interval." ! ! Alexander Godunov ! Companion website, 2026 !======================================================================= program main implicit none double precision :: g double precision :: a, b, integral, exact integer :: n, i integer, dimension(4) :: nvalues external :: g a = 0.0d0 b = 1.0d0 exact = 0.5d0 nvalues = (/ 8, 16, 32, 64 /) write(*,'(a)') 'Improper integral: infinite interval' write(*,'(a)') 'Integral of exp(-x) sin(x) from 0 to infinity' write(*,'(a,es15.7)') 'Exact value = ', exact write(*,*) write(*,'(a)') ' n Gauss-Legendre' write(*,'(a)') '--------------------------------' do i = 1, size(nvalues) n = nvalues(i) call gauss_legendre(g, a, b, n, integral) write(*,'(i9,es21.7)') n, integral end do end program main double precision function f(x) !--------------------------------------------------------------------- ! Original function to be integrated on the semi-infinite interval. ! ! infinity ! / ! | f(x) dx ! / ! 0 ! ! To use this program for another problem, replace only the line ! defining f(x) below. The transformation and quadrature routines ! do not need to be changed. ! ! Input: ! x - original integration variable, 0 <= x < infinity ! ! Output: ! f - value of the original integrand !--------------------------------------------------------------------- implicit none double precision, intent(in) :: x f = exp(-x)*sin(x) end function f double precision function g(t) !--------------------------------------------------------------------- ! Transformed integrand on the finite interval 0 <= t < 1. ! ! The semi-infinite interval is mapped to [0,1) using ! ! x = t/(1-t), dx = dt/(1-t)^2 . ! ! Therefore ! ! g(t) = f(x)/(1-t)^2, ! ! where x = t/(1-t). ! ! This function performs only the variable transformation. The ! original integrand is defined separately in function f(x), so a ! user can change the physical or mathematical problem by modifying ! f(x) without changing this routine. ! ! Input: ! t - transformed integration variable, 0 <= t < 1 ! ! Output: ! g - transformed integrand used by the finite-interval quadrature ! ! The point t = 1 corresponds to x = infinity. Gauss-Legendre ! quadrature does not sample the endpoints, so f is never evaluated ! at x = infinity. !--------------------------------------------------------------------- implicit none double precision :: f double precision, intent(in) :: t double precision :: x, one_minus_t external :: f one_minus_t = 1.0d0 - t x = t/one_minus_t g = f(x)/(one_minus_t*one_minus_t) end function g subroutine gauss_legendre(f, a, b, n, integral) !--------------------------------------------------------------------- ! n-point Gauss-Legendre quadrature for integration of f(x) on [a,b]. ! ! Input: ! f - function to integrate ! a - lower integration limit ! b - upper integration limit ! n - number of Gauss-Legendre points ! ! Output: ! integral - numerical approximation to the integral ! ! Method: ! The nodes are the roots of the Legendre polynomial P_n(x). ! They are found by Newton iteration. The corresponding weights ! are calculated from P'_n(x). Symmetry is used so that only half ! of the roots need to be computed explicitly. ! ! The nodes and weights are defined on [-1,1] and are mapped to ! the interval [a,b]. ! ! Based on the original Gauss quadrature routines written by ! Alexander Godunov, October 2009. !--------------------------------------------------------------------- implicit none double precision :: f double precision, intent(in) :: a, b integer, intent(in) :: n double precision, intent(out) :: integral external :: f double precision, parameter :: pi = 3.14159265358979323846d0 double precision, parameter :: tol = 1.0d-14 integer, parameter :: max_iter = 100 double precision :: z, zold double precision :: p0, p1, p2, dp double precision :: weight double precision :: midpoint, halfwidth double precision :: xleft, xright integer :: i, j, m, iter if (n < 1) then write(*,'(a)') 'gauss_legendre: n must be positive.' stop end if midpoint = (a+b)/2.0d0 halfwidth = (b-a)/2.0d0 integral = 0.0d0 ! Because the roots are symmetric, only half need to be calculated. m = (n+1)/2 do i = 1, m ! Initial approximation to a root of P_n. z = cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0)) ! Newton iteration. do iter = 1, max_iter p0 = 1.0d0 p1 = z if (n == 1) then p2 = p1 else do j = 2, n p2 = ((2.0d0*dble(j)-1.0d0)*z*p1 & - (dble(j)-1.0d0)*p0)/dble(j) p0 = p1 p1 = p2 end do end if if (n == 1) then dp = 1.0d0 else dp = dble(n)*(z*p1-p0)/(z*z-1.0d0) end if zold = z z = zold - p1/dp if (abs(z-zold) <= tol) exit end do if (iter > max_iter) then write(*,'(a)') 'gauss_legendre: Newton iteration did not converge.' stop end if ! Re-evaluate P_n and P_{n-1} at the converged root. p0 = 1.0d0 p1 = z if (n == 1) then dp = 1.0d0 else do j = 2, n p2 = ((2.0d0*dble(j)-1.0d0)*z*p1 & - (dble(j)-1.0d0)*p0)/dble(j) p0 = p1 p1 = p2 end do dp = dble(n)*(z*p1-p0)/(z*z-1.0d0) end if weight = 2.0d0/((1.0d0-z*z)*dp*dp) xleft = midpoint - halfwidth*z xright = midpoint + halfwidth*z if (abs(z) <= tol) then integral = integral + weight*f(midpoint) else integral = integral + weight*(f(xleft) + f(xright)) end if end do integral = halfwidth*integral end subroutine gauss_legendre