gaussjacobiquad Module Reference

Overall driver for Gauss-Jacobi quadrature. More...

Functions/Subroutines
subroutine	gauss_jacobi (npts, alpha, beta, x, wts, method)
	Compute the Gauss-Jacobi quadrature nodes and weights. More...

Detailed Description

Overall driver for Gauss-Jacobi quadrature.

Function/Subroutine Documentation

gauss_jacobi()

subroutine gaussjacobiquad::gauss_jacobi	(	integer, intent(in)	npts,
		real(dp), intent(in)	alpha,
		real(dp), intent(in)	beta,
		real(dp), dimension(npts), intent(out)	x,
		real(dp), dimension(npts), intent(out)	wts,
		character(len=:), intent(in), allocatable	method
	)

Compute the Gauss-Jacobi quadrature nodes and weights.

This subroutine calculates the Gauss-Jacobi quadrature nodes and weights for the given parameters \(\alpha\) and \(\beta\), using the specified method. Gauss-Jacobi quadrature is used to approximate integrals of the form: [ \int_{-1}^{1} (1 - x)^\alpha (1 + x)^\beta f(x) \,dx \approx \sum_{i=1}^{npts} wts_i f(x_i) ] where the weights and nodes are calculated with the Jacobi polynomial, which is defined as: [ P_n^{(\alpha, \beta)}(x) = \frac{(\alpha + 1)_n}{n!} \sum_{k=0}^n \binom{n}{k} \frac{(\beta + 1)_{n-k}}{(n-k)!} \left( \frac{x-1}{2} \right)^k \left( \frac{x+1}{2} \right)^{n-k} ]

Parameters

[in]	npts	Number of quadrature points.
[in]	alpha	Parameter alpha in the Jacobi polynomial. Must be greater than -1.
[in]	beta	Parameter beta in the Jacobi polynomial. Must be greater than -1.
[out]	x	Quadrature nodes.
[out]	wts	Quadrature weights.
[in]	method	Method used for calculation. Supported methods are "rec" and "gw".

Definition at line 50 of file GaussJacobiQuad.f90.

    integer, intent(in) :: npts
    real(dp), intent(in) :: alpha, beta
    real(dp), intent(out) :: x(npts), wts(npts)
    character(len=:), allocatable, intent(in) :: method
 
    if (npts <= 0) then
        error stop "Number of points must be positive"
    end if
 
    if (alpha <= -1.0_dp) then
        error stop "alpha must be greater than -1"
    end if
 
    if (beta <= -1.0_dp) then
 
        error stop "beta must be greater than -1"
    end if
 
    select case (trim(method))
    case ("rec") ! Fails at high beta
        call gauss_jacobi_rec(npts, alpha, beta, x, wts)
    case ("gw") ! Accurate for high beta
        call gauss_jacobi_gw(npts, alpha, beta, x, wts)
    case ("algo665")
        call gauss_jacobi_algo665(npts, alpha, beta, x, wts)
    case default
        print*,"Error: Unknown method specified:", method
        print*,"Supported methods: 'rec', 'gw', 'algo665'"
        error stop
    end select