module gauss_quadrature_m use constants_m, only: PI use precision_m, only: DP, FP implicit none private public :: gauss_legendre, gauss_laguerre contains pure subroutine gauss_legendre(x, w) !! Calculates the nodes and weights to perform Gauss-Legendre !! quadrature on the interval [-1, 1] !! !! NOTE: This subroutine has been tested for array lengths up to 64 !! and delivers an accuracy of 1e-14. real(kind = FP), intent(out) :: x(:) !! Quadrature nodes real(kind = FP), intent(out) :: w(:) !! Quadrature weights real(kind = DP), dimension(:), allocatable :: x_dp, w_dp associate(n => size(x)) allocate(x_dp(n)) allocate(w_dp(n)) select case (n) case (1) x_dp = 0 w_dp = 2 case (2) x_dp(1) = -sqrt(1.0_DP/3.0_DP) x_dp(2) = -x_dp(1) w_dp = 1 case default block ! The following implementation has been adapted from Numerical ! Recipies in Fortran 77 integer :: i, j, m real(kind = DP) :: p1, p2, p3, pp, z, z1 real(kind = DP), parameter :: tolerance = 1.0e-15_DP m = (n+1)/2 ! The roots are symmetric in the interval, so we only have to ! find half of them do i=1,m ! Initial guess z = cos(PI * (i - 0.25_DP) / (n + 0.5_DP)) z1 = 0.0 ! Main loop of refinement by Newton’s method do while(abs(z - z1) > tolerance) p1 = 1.0_DP p2 = 0.0_DP do j = 1, n ! Loop up the recurrence relation to get the ! Legendre polynomial evaluated at z p3 = p2 p2 = p1 p1 = ((2.0_DP * j - 1.0_DP) * z * p2 & - (j - 1.0_DP) * p3) / j end do ! p1 is now the desired Legendre polynomial. We next ! compute pp, its derivative, by a standard relation ! involving also p2, the polynomial of one lower order. pp = n * (z * p1 - p2) / (z * z - 1.0_DP) z1 = z z = z1 - p1 / pp ! Newton’s method end do x_dp(i) = -z x_dp(n+1-i) = z w_dp(i) = 2.0_DP / ((1.0_DP - z * z) * pp * pp) w_dp(n+1-i) = w_dp(i) end do end block end select end associate x = real(x_dp, kind = FP) w = real(w_dp, kind = FP) deallocate(x_dp, w_dp) end subroutine subroutine gauss_laguerre(x, w) !! Calculates the nodes and weights to perform Gauss-Laguerre !! quadrature on the interval [0, infinity) !! !! NOTE: This subroutine has been tested for array lengths up to 64 !! and delivers an accuracy of 1e-12 for 64 array elements in !! double precision. For fewer elements the accuracy is !! improved. real(kind = FP), intent(out) :: x(:) !! Quadrature nodes real(kind = FP), intent(out) :: w(:) !! Quadrature weights real(kind = DP), dimension(:), allocatable :: x_dp, w_dp associate(n => size(x)) allocate(x_dp(n)) allocate(w_dp(n)) select case (n) case (1) x_dp = 1.0_DP w_dp = 1.0_DP case (2) x_dp(1) = 2.0_DP - sqrt(2.0_DP) x_dp(2) = 2.0_DP + sqrt(2.0_DP) w_dp(1) = 0.25_DP * (2.0_DP + sqrt(2.0_DP)) w_dp(2) = 0.25_DP * (2.0_DP - sqrt(2.0_DP)) case default block ! The following implementation has been adapted from GENE integer :: nr, i, k, j, it integer, parameter :: max_it_count = 40 real(kind = DP) :: hn, pf, pd, z, z0, p, f0, f1, fd, q, wp, gd real(kind = DP), parameter :: tolerance = 1.0e-15_DP hn = 1.0_DP / n pf = 0.0_DP pd = 0.0_DP do nr = 1, n ! Initial guess z = hn if (nr > 1) z = x_dp(nr - 1) + hn * nr**1.27_DP ! Start the Newton iteration it = 0 10 it = it + 1 z0 = z p = 1.0_DP do i = 1, nr - 1 p = p * (z - x_dp(i)) enddo f0 = 1.0_DP f1 = 1.0_DP - z ! Loop up the recurrence relation to get the Laguerre ! polynomial evaluated at z do k = 2, n pf = ((2.0_DP * k - 1.0_DP - z) * f1 & - (k - 1.0_DP) * f0) / k pd = k / z * (pf - f1) f0 = f1 f1 = pf enddo ! pf is the desired Laguerre polynomial and pd its ! derivative fd = pf / p q = 0.0_DP do i = 1, nr - 1 wp = 1.0_DP do j= 1, nr - 1 if (j /= i) wp = wp * (z - x_dp(j)) enddo q = q + wp enddo gd = (pd - q * fd) / p ! Newton's method z = z - fd / gd if (it <= max_it_count & .and. abs((z - z0) / z) > tolerance) go to 10 x_dp(nr) = z w_dp(nr) = 1.0_DP / (z * pd * pd) enddo end block end select end associate x = real(x_dp, kind = FP) w = real(w_dp, kind = FP) deallocate(x_dp, w_dp) end subroutine end module