get_intrvl_weights Subroutine

private subroutine get_intrvl_weights(llim, ulim, nodes, wgts)

Used by get_weights to find the Lagrange quadrature weights for the given grid.

Suppose we wish to calculate Int(f[x] dx) over the range [llim, ulim]. We could approximate f[x] by the Lagrange interpolating polynomial --> f[x] ~ Sum_i{ f[x_i] l_i(x)} where. the Lagrange polynomial is given by l_i(x) = Pi_j{ (x - x_i)/(x_i - x_j) }. Hence our integral can be written as I = Int(f[x] dx) ~ Sum_i{ f[x_i] Int(Pi_j{ (x - x_i)/(x_i - x_j) })} so I ~ Sum_i{f[x_i] c_i} with c_i = Int(Pi_j{ (x - x_i)/(x_i - x_j) }). Here we are calculating the c_i values which are passed out as wgts corresponding to the set of positions {x_i} set by nodes. To evaluate the integral of the Lagrange polynomial coefficients we adopt Gauss-Legendre quadrature, which allows highly accurate evaluation of polynomials of degree 2*N-1. Given the need to find M weights we find we only need Gauss-Legendre of order N = (M + 1)/2. Here we actually use degree (M/2) + 1.

One limitation of Lagrange polynomials is their lack of numerical stability at high order (i.e. when M is large, typically > 20 is quoted as a limit). The routine get_weights allows the integration domain to be sub-divided in order to limit the Lagrange order and reduce the variation in weights. This is partially controlled by the input nmax.

Arguments

  subroutine get_intrvl_weights (llim, ulim, nodes, wgts)
    use gauss_quad, only: get_legendre_grids_from_cheb
    
    implicit none
    
    ! llim (ulim) is lower (upper) limit of integration
    real, intent (in) :: llim, ulim
    real, dimension (:), intent (in) :: nodes
    real, dimension (:), intent (in out) :: wgts
    
    ! stuff needed to do guassian quadrature 
    real, dimension (:), allocatable :: gnodes, gwgts, omprod
    integer :: ix, iw

    allocate (gnodes(size(nodes)/2+1), gwgts(size(wgts)/2+1), omprod(size(nodes)/2+1))
    
    call get_legendre_grids_from_cheb (llim, ulim, gnodes, gwgts)

    do iw=1,size(wgts)
       omprod = 1.0
       
       do ix=1,size(nodes)
          if (ix /= iw) omprod = omprod*(gnodes - nodes(ix))/(nodes(iw) - nodes(ix))
       end do
       
       do ix=1,size(gwgts)
          wgts(iw) = wgts(iw) + omprod(ix)*gwgts(ix)
       end do
    end do

    deallocate (gnodes, gwgts, omprod)
       
  end subroutine get_intrvl_weights