! gaussquad.f90: Module for gaussian quadrature rules ! http://infty.us/gaussquad/gaussquad.html ! v1.0 ! ! Copyright (c) 2012 Christopher N. Gilbreth ! ! Permission is hereby granted, free of charge, to any person obtaining a copy ! of this software and associated documentation files (the "Software"), to deal ! in the Software without restriction, including without limitation the rights ! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell ! copies of the Software, and to permit persons to whom the Software is ! furnished to do so, subject to the following conditions: ! ! The above copyright notice and this permission notice shall be included in all ! copies or substantial portions of the Software. ! ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE ! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, ! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE ! SOFTWARE. module gaussquad implicit none ! References: ! 1. G. H. Golub and J. H. Welsch, Math. Comp. 23, 221 (1969) ! 2. Digital Library of Mathematical Functions, Ch. 3.5, 18.9 integer, parameter :: rk = selected_real_kind(p=15) contains subroutine cpquad(n,alpha,name,w,x) ! Compute the weights w_i and abscissas x_i for an n-point Gaussian ! quadrature rule for classical polynomials of type "name". ! Input: ! n: Degree of quadrature (= number of points) ! alpha: Parameter associated with the polynomials ! name: Name of polynomials: "Legendre", "Laguerre", "Hermite" ! Output: ! w: Weights of the quadrature rule ! x: Abscissas of the quadrature rule. ! ! Notes: ! The weights w(i) and abscissas x(i) are computed so that ! ∫W(x) f(x) dx = ∑_i w(i) f(x(i)) ! exactly for polynomials f of degree ≤ 2n-1. ! ! Here W(x) is the wieght function associated with the orthogonal ! polynomials p_i(x), i.e. ! ∫W(x) p_i(x) p_j(x) dx = C_i δ_{i,j}. ! Each type of polynomial corresponds to a particlar weight function W. ! ! name weight function W(x) ! ---- -------------------- ! Legendre 1 ! Laguerre x^{alpha} exp(-x) ! Hermite exp(-x^2) ! ! The parameter alpha only affects the Laguerre calculation. implicit none integer, intent(in) :: n real(rk), intent(in) :: alpha character(len=*), intent(in) :: name real(rk), intent(out) :: w(n), x(n) real(rk) :: diag(n), subdiag(n-1), work(max(1,2*n-2)), mu0 integer :: info, i real(rk), allocatable :: P(:,:) allocate(P(n,n)) ! Compute the symmetric tridiagonal matrix "J" for the given type of polynomial call cpj(n,alpha,name,diag,subdiag) ! Diagonalize (using LAPACK) call DSTEV('V', n, diag, subdiag, P, n, work, info) if (info .ne. 0) stop "Error diagonalizing in cpquad" ! Eigenvalues are the abscissas, i.e. the zeroes of p_n(x) x = diag ! Calculate the zeroth moment of the weight function W(x) call cpzm(alpha,name,mu0) ! Weights (Eq. 2.6 in Golub & Welsch) do i=1,n w(i) = P(1,i)**2 * mu0 end do end subroutine cpquad subroutine cpzm(alpha,name,mu0) ! Calculate ! μ_0 = ∫W(x)dx, ! the zeroth moment of the wieght function for a given type of classical ! orthgonal polynomial. ! Input: ! alpha: Parameter associated with the polynomials ! name: Name of the polynomials ! Output: ! mu0: The zeroth moment, as above. implicit none real(rk), intent(in) :: alpha character(len=*), intent(in) :: name real(rk), intent(out) :: mu0 if (name == "Legendre") then mu0 = 2 else if (name == "Laguerre") then mu0 = gamma(alpha + 1) else if (name == "Hermite") then mu0 = sqrt(2*asin(1._rk)) else stop "Unknown polynomial in cpzm" end if end subroutine cpzm subroutine cpj(n,alpha,name,diag,subdiag) ! Calculate the symmetric tridiagonal matrix J associated with the classical ! orthogonal polynimals given by 'name'. ! Input: ! n: Degree (= number of roots) ! alpha: Parameter associated with the polynomials ! name: Name of polynimials: "legendre", "laguerre", "hermite" ! Output: ! diag: Diagonal matrix elements of J ! subdiag: Subdiagonal matrix elements of J ! Notes: ! Reference: G. H. Golub and J. H. Welsch, Math. Comp. 23, 221 (1969) implicit none integer, intent(in) :: n real(rk), intent(in) :: alpha character(len=*), intent(in) :: name real(rk), intent(out) :: diag(n), subdiag(n-1) integer :: i ! In general, if ! P_n(x) = (a_n x + b_n) P_{n-1}(x) - c_n P_{n-2}(x) ! Then the n diagonal elements of the matrix are ! α_i = -b_i/a_i, i=1,...,n ! And the n-1 subdiagonal elements are ! β_i = Sqrt(c_{i+1}/(a_i a_{i+1})), i=1,...,n-1. if (name == "Legendre") then ! Legendre polynomials: ! P_n(x) = (a_n x + b_n) P_{n-1}(x) - c_n P_{n-2}(x) ! where ! a_n = (2n-1)/n ! b_n = 0 ! c_n = (n-1)/n. ! So ! α_i = -b_i/a_i = 0 ! β_i = Sqrt(i^2/[(2i-1)(2i+1)]) diag = 0._rk do i=1,n-1 subdiag(i) = sqrt(i**2 / ((2*i-1._rk)*(2*i+1._rk))) end do else if (name == "Laguerre") then ! Associated Laguerre polynomials: ! L^{α}_n(x) = (a_n x + b_n) L^{α}_{n-1}(x) - c_n L^{α}_{n-2}(x) ! where ! a_n = -1/n ! b_n = (α + 2n-1)/n ! c_n = (α + n-1)/n ! So ! α_i = (α + 2i-1) ! β_i = sqrt((α+i)i) do i=1,n diag(i) = alpha+2*i-1 end do do i=1,n-1 subdiag(i) = sqrt((alpha+i)*i) end do else if (name == "Hermite") then ! Hermite polynomials: ! H_n(x) = (a_n x + b_n) H_{n-1}(x) - c_n H_{n-2}(x) ! where ! a_n = 2 ! b_n = 0 ! c_n = 2*(n-1) ! So ! α_i = 0 ! β_i = sqrt(i/2) diag = 0._rk do i=1,n-1 subdiag(i) = sqrt(i/2._rk) end do else stop "Unknown polynomial in cpj" end if end subroutine cpj end module gaussquad