gaussquad (version 1.0-1)

jacobi.p.quadrature: Perform Gauss Jacobi quadrature

Description

This function evaluates the integral of the given function between the lower and upper limits using the weight and abscissa values specified in the rule data frame. The quadrature formula uses the weight function for Jacobi P polynomials.

Usage

jacobi.p.quadrature(functn, rule, alpha = 0, beta = 0, lower = -1, upper = 1, 
weighted = TRUE,  ...)

Arguments

functn

an R function which should take a numeric argument x and possibly some parameters. The function returns a numerical vector value for the given argument x.

rule

a data frame containing the order $n$ ultraspherical quadrature rule

alpha

numeric value for the first Jacobi polynomial parameter

beta

numeric value for the second Jacobi polynomial parameter

lower

numeric value for the lower limit of the integral with a default value of -1

upper

numeric value for the upper limit of the integral with a default value of 1

weighted

a boolean value which if true causes the ultraspherical weight function to be included in the integrand

...

other arguments passed to the give function

Value

• The value of definite integral evaluated using Gauss Jacobi quadrature

Details

The rule argument corresponds to an order $n$ Jacobi polynomial, weight function and interval $\left[ { - 1,1} \right]$. The lower and upper limits of the integral must be finite.

References

Abramowitz, M. and I. A. Stegun, 1968. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover Publications, Inc., New York. Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, 1992. Numerical Recipes in C, Cambridge University Press, Cambridge, U.K. Stroud, A. H., and D. Secrest, 1966. Gaussian Quadrature Formulas, Prentice-Hall, Englewood Cliffs, NJ.

See Also

jacobi.p.quadrature.rules

Examples

###
### this example evaluates the quadrature function for
### the Jacobi P polynomials.  it computes the integral
### of the product for all pairs of orthogonal polynomials
### from order 0 to order 5.  the results are compared to
### the diagonal matrix of the inner products for the
### polynomials
###
### construct a list of the Jacobi P polynomials
### of orders 0 to 5
### alpha is 2 and beta is 2
###
alpha <- 2
beta <- 2
p.list <- jacobi.p.polynomials( 5, alpha, beta )
###
### function to construct the polynomial products by column
###
by.column.products <- function( c, p.list, p.p.list )
{
###
### function to construct the polynomial products by row
###
    by.row.products <- function( r, c, p.list )
    {
        row.column.product <- p.list[[r]] * p.list[[c]]
        return (row.column.product )
    }
    row.list <- lapply( 1:6, by.row.products, c, p.list )
    return( row.list )
}
###
### construct the two dimensional list of pair products
### of polynomials
###
p.p.products <- lapply( 1:6, by.column.products, p.list )
###
### function construct the polynomial functions by column
###
by.column.functions <- function( c, p.p.products )
{
###
### function to construct the polynomial functions by row
###
    by.row.functions <- function( r, c, p.p.products )
    {
        row.column.function <- as.function( p.p.products[[r]][[c]] )
        return( row.column.function )
    }
    row.list <- lapply( 1:6, by.row.functions, c, p.p.products )
    return( row.list )
}
###
### compute the two dimensional list of functions
### corresponding to the polynomial products in
### the two dimensional list p.p.products
###
p.p.functions <- lapply( 1:6, by.column.functions, p.p.products )
###
### get the rule table for the order 6 polynomial
###
rules <- jacobi.p.quadrature.rules( 6, alpha, beta )
order.6.rule <- rules[[6]]
###
### function to compute the integral of the polynomials by column
###
by.column.integrals <- function( c, p.p.functions )
{
###
### function to compute the integral of the polynomials by row
###
    by.row.integrals <- function( r, c, p.p.functions )
    {
        row.column.integral <- jacobi.p.quadrature(
            p.p.functions[[r]][[c]], order.6.rule, alpha, beta )
        return( row.column.integral )
    }
    row.vector <- sapply( 1:6, by.row.integrals, c, p.p.functions )
    return( row.vector )
}
###
### construct the square symmetric matrix containing
### the definite integrals over the default limits
### corresponding to the two dimensional list of
### polynomial functions
###
p.p.integrals <- sapply( 1:6, by.column.integrals, p.p.functions )
###
### construct the diagonal matrix with the inner products
### of the orthogonal polynomials on the diagonal
###
p.p.inner.products <- diag( jacobi.p.inner.products( 5, alpha, beta ) )
print( p.p.integrals )
print( p.p.inner.products )