Functions
Polylib Namespace Reference

The namespace associated with the the Polylib library (Polylib introduction) More...

Functions

static void Jacobz (const int n, double *z, const double alpha, const double beta)
 Calculate the n zeros, z, of the Jacobi polynomial, i.e. More...
 
static void TriQL (const int n, double *d, double *e, double **z)
 QL algorithm for symmetric tridiagonal matrix. More...
 
double gammaF (const double x)
 Calculate the Gamma function , $ \Gamma(n)$, for integer. More...
 
static void RecCoeff (const int n, double *a, double *b, const double alpha, const double beta)
 The routine finds the recurrence coefficients a and. More...
 
void JKMatrix (int n, double *a, double *b)
 Calcualtes the Jacobi-kronrod matrix by determining the. More...
 
void chri1 (int n, double *a, double *b, double *a0, double *b0, double z)
 Given a weight function $w(t)$ through the first n+1. More...
 
void zwgj (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Jacobi zeros and weights. More...
 
void zwgrjm (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Radau-Jacobi zeros and weights with end point at z=-1. More...
 
void zwgrjp (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Radau-Jacobi zeros and weights with end point at z=1. More...
 
void zwglj (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1. More...
 
void zwgk (double *z, double *w, const int npt, const double alpha, const double beta)
 Gauss-Kronrod-Jacobi zeros and weights. More...
 
void zwrk (double *z, double *w, const int npt, const double alpha, const double beta)
 Gauss-Radau-Kronrod-Jacobi zeros and weights. More...
 
void zwlk (double *z, double *w, const int npt, const double alpha, const double beta)
 Gauss-Lobatto-Kronrod-Jacobi zeros and weights. More...
 
void Dgj (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix and its transpose associated. More...
 
void Dgrjm (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix and its transpose associated. More...
 
void Dgrjp (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix associated with the. More...
 
void Dglj (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix associated with the. More...
 
double hgj (const int i, const double z, const double *zgj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through. More...
 
double hgrjm (const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through the. More...
 
double hgrjp (const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through the. More...
 
double hglj (const int i, const double z, const double *zglj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through the. More...
 
void Imgj (double *im, const double *zgj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Jacobi points to an. More...
 
void Imgrjm (double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Radau-Jacobi points. More...
 
void Imgrjp (double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Radau-Jacobi points. More...
 
void Imglj (double *im, const double *zglj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Lobatto-Jacobi points. More...
 
void jacobfd (const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
 Routine to calculate Jacobi polynomials, $ P^{\alpha,\beta}_n(z) $, and their first derivative, $ \frac{d}{dz} P^{\alpha,\beta}_n(z) $. More...
 
void jacobd (const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
 Calculate the derivative of Jacobi polynomials. More...
 
void JacZeros (const int n, double *a, double *b, const double alpha, const double beta)
 Zero and Weight determination through the eigenvalues and eigenvectors of a tridiagonal. More...
 

Detailed Description

The namespace associated with the the Polylib library (Polylib introduction)

Function Documentation

void Polylib::chri1 ( int  n,
double *  a,
double *  b,
double *  a0,
double *  b0,
double  z 
)

Given a weight function $w(t)$ through the first n+1.

coefficients a and b of its orthogonal polynomials

this routine generates the first n recurrence coefficients for the orthogonal

polynomials relative to the modified weight function $(t-z)w(t)$.

The result will be placed in the array a0 and b0.

Here is the caller graph for this function:

LIB_UTILITIES_EXPORT void Polylib::Dgj ( double *  D,
const double *  z,
const int  np,
const double  alpha,
const double  beta 
)

Compute the Derivative Matrix and its transpose associated.

with the Gauss-Jacobi zeros.

  • Compute the derivative matrix, d, and its transpose, dt,

associated with the n_th order Lagrangian interpolants through the

np Gauss-Jacobi points z such that
$ \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) $

Here is the call graph for this function:

LIB_UTILITIES_EXPORT void Polylib::Dglj ( double *  D,
const double *  z,
const int  np,
const double  alpha,
const double  beta 
)

Compute the Derivative Matrix associated with the.

Gauss-Lobatto-Jacobi zeros.

  • Compute the derivative matrix, d, associated with the n_th

order Lagrange interpolants through the np

Gauss-Lobatto-Jacobi points z such that
$ \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) $

Here is the call graph for this function:

LIB_UTILITIES_EXPORT void Polylib::Dgrjm ( double *  D,
const double *  z,
const int  np,
const double  alpha,
const double  beta 
)

Compute the Derivative Matrix and its transpose associated.

with the Gauss-Radau-Jacobi zeros with a zero at z=-1.

  • Compute the derivative matrix, d, associated with the n_th

order Lagrangian interpolants through the np Gauss-Radau-Jacobi

points z such that
$ \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) $

Here is the call graph for this function:

LIB_UTILITIES_EXPORT void Polylib::Dgrjp ( double *  D,
const double *  z,
const int  np,
const double  alpha,
const double  beta 
)

Compute the Derivative Matrix associated with the.

Gauss-Radau-Jacobi zeros with a zero at z=1.

  • Compute the derivative matrix, d, associated with the n_th

order Lagrangian interpolants through the np Gauss-Radau-Jacobi

points z such that
$ \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) $

Here is the call graph for this function:

double Polylib::gammaF ( const double  x)

Calculate the Gamma function , $ \Gamma(n)$, for integer.

values and halves.

Determine the value of $\Gamma(n)$ using:

$ \Gamma(n) = (n-1)! \mbox{ or } \Gamma(n+1/2) = (n-1/2)\Gamma(n-1/2)$

where $ \Gamma(1/2) = \sqrt(\pi)$

Here is the caller graph for this function:

LIB_UTILITIES_EXPORT double Polylib::hgj ( const int  i,
const double  z,
const double *  zgj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through.

the np Gauss-Jacobi points zgj at the arbitrary location z.

  • $ -1 \leq z \leq 1 $
  • Uses the defintion of the Lagrangian interpolant:
    %

$ \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{P_{np}^{\alpha,\beta}(z)} {[P_{np}^{\alpha,\beta}(z_j)]^\prime (z-z_j)} & \mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} $

Here is the call graph for this function:

Here is the caller graph for this function:

LIB_UTILITIES_EXPORT double Polylib::hglj ( const int  i,
const double  z,
const double *  zglj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through the.

np Gauss-Lobatto-Jacobi points zgrj at the arbitrary location

z.

  • $ -1 \leq z \leq 1 $
  • Uses the defintion of the Lagrangian interpolant:
    %

$ \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{(1-z^2) P_{np-2}^{\alpha+1,\beta+1}(z)} {((1-z^2_j) [P_{np-2}^{\alpha+1,\beta+1}(z_j)]^\prime - 2 z_j P_{np-2}^{\alpha+1,\beta+1}(z_j) ) (z-z_j)}&\mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} $

Here is the call graph for this function:

Here is the caller graph for this function:

LIB_UTILITIES_EXPORT double Polylib::hgrjm ( const int  i,
const double  z,
const double *  zgrj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through the.

np Gauss-Radau-Jacobi points zgrj at the arbitrary location

z. This routine assumes zgrj includes the point -1.

  • $ -1 \leq z \leq 1 $
  • Uses the defintion of the Lagrangian interpolant:
    %

$ \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{(1+z) P_{np-1}^{\alpha,\beta+1}(z)} {((1+z_j) [P_{np-1}^{\alpha,\beta+1}(z_j)]^\prime + P_{np-1}^{\alpha,\beta+1}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} $

Here is the call graph for this function:

Here is the caller graph for this function:

LIB_UTILITIES_EXPORT double Polylib::hgrjp ( const int  i,
const double  z,
const double *  zgrj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through the.

np Gauss-Radau-Jacobi points zgrj at the arbitrary location

z. This routine assumes zgrj includes the point +1.

  • $ -1 \leq z \leq 1 $
  • Uses the defintion of the Lagrangian interpolant:
    %

$ \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{(1-z) P_{np-1}^{\alpha+1,\beta}(z)} {((1-z_j) [P_{np-1}^{\alpha+1,\beta}(z_j)]^\prime - P_{np-1}^{\alpha+1,\beta}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} $

Here is the call graph for this function:

Here is the caller graph for this function:

LIB_UTILITIES_EXPORT void Polylib::Imgj ( double *  im,
const double *  zgj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Jacobi points to an.

arbitrary distribution at points zm

  • Computes the one-dimensional interpolation matrix, im, to

interpolate a function from at Gauss-Jacobi distribution of nz

zeros zgrj to an arbitrary distribution of mz points zm, i.e.
$ u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) $

Here is the call graph for this function:

LIB_UTILITIES_EXPORT void Polylib::Imglj ( double *  im,
const double *  zglj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Lobatto-Jacobi points.

to an arbitrary distrubtion at points zm

  • Computes the one-dimensional interpolation matrix, im, to

interpolate a function from at Gauss-Lobatto-Jacobi distribution of

nz zeros zgrj (where zgrj[0]=-1) to an arbitrary

distribution of mz points zm, i.e.


$ u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) $

Here is the call graph for this function:

LIB_UTILITIES_EXPORT void Polylib::Imgrjm ( double *  im,
const double *  zgrj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Radau-Jacobi points.

(including z=-1) to an arbitrary distrubtion at points zm

  • Computes the one-dimensional interpolation matrix, im, to

interpolate a function from at Gauss-Radau-Jacobi distribution of

nz zeros zgrj (where zgrj[0]=-1) to an arbitrary

distribution of mz points zm, i.e.


$ u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) $

Here is the call graph for this function:

LIB_UTILITIES_EXPORT void Polylib::Imgrjp ( double *  im,
const double *  zgrj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Radau-Jacobi points.

(including z=1) to an arbitrary distrubtion at points zm

  • Computes the one-dimensional interpolation matrix, im, to

interpolate a function from at Gauss-Radau-Jacobi distribution of

nz zeros zgrj (where zgrj[nz-1]=1) to an arbitrary

distribution of mz points zm, i.e.


$ u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) $

Here is the call graph for this function:

LIB_UTILITIES_EXPORT void Polylib::jacobd ( const int  np,
const double *  z,
double *  polyd,
const int  n,
const double  alpha,
const double  beta 
)

Calculate the derivative of Jacobi polynomials.

  • Generates a vector poly of values of the derivative of the

n th order Jacobi polynomial $ P^(\alpha,\beta)_n(z)$ at the

np points z.

  • To do this we have used the relation


$ \frac{d}{dz} P^{\alpha,\beta}_n(z) = \frac{1}{2} (\alpha + \beta + n + 1) P^{\alpha,\beta}_n(z) $

  • This formulation is valid for $ -1 \leq z \leq 1 $

Here is the call graph for this function:

Here is the caller graph for this function:

LIB_UTILITIES_EXPORT void Polylib::jacobfd ( const int  np,
const double *  z,
double *  poly_in,
double *  polyd,
const int  n,
const double  alpha,
const double  beta 
)

Routine to calculate Jacobi polynomials, $ P^{\alpha,\beta}_n(z) $, and their first derivative, $ \frac{d}{dz} P^{\alpha,\beta}_n(z) $.

  • This function returns the vectors poly_in and poly_d

containing the value of the $ n^th $ order Jacobi polynomial

$ P^{\alpha,\beta}_n(z) \alpha > -1, \beta > -1 $ and its

derivative at the np points in z[i]

  • If poly_in = NULL then only calculate derivatice
  • If polyd = NULL then only calculate polynomial
  • To calculate the polynomial this routine uses the recursion

relationship (see appendix A ref [4]) :

$ \begin{array}{rcl} P^{\alpha,\beta}_0(z) &=& 1 \\ P^{\alpha,\beta}_1(z) &=& \frac{1}{2} [ \alpha-\beta+(\alpha+\beta+2)z] \\ a^1_n P^{\alpha,\beta}_{n+1}(z) &=& (a^2_n + a^3_n z) P^{\alpha,\beta}_n(z) - a^4_n P^{\alpha,\beta}_{n-1}(z) \\ a^1_n &=& 2(n+1)(n+\alpha + \beta + 1)(2n + \alpha + \beta) \\ a^2_n &=& (2n + \alpha + \beta + 1)(\alpha^2 - \beta^2) \\ a^3_n &=& (2n + \alpha + \beta)(2n + \alpha + \beta + 1) (2n + \alpha + \beta + 2) \\ a^4_n &=& 2(n+\alpha)(n+\beta)(2n + \alpha + \beta + 2) \end{array} $

  • To calculate the derivative of the polynomial this routine uses

the relationship (see appendix A ref [4]) :

$ \begin{array}{rcl} b^1_n(z)\frac{d}{dz} P^{\alpha,\beta}_n(z)&=&b^2_n(z)P^{\alpha,\beta}_n(z) + b^3_n(z) P^{\alpha,\beta}_{n-1}(z) \hspace{2.2cm} \\ b^1_n(z) &=& (2n+\alpha + \beta)(1-z^2) \\ b^2_n(z) &=& n[\alpha - \beta - (2n+\alpha + \beta)z]\\ b^3_n(z) &=& 2(n+\alpha)(n+\beta) \end{array} $

  • Note the derivative from this routine is only valid for -1 < z < 1.

Here is the caller graph for this function:

static void Polylib::Jacobz ( const int  n,
double *  z,
const double  alpha,
const double  beta 
)
static

Calculate the n zeros, z, of the Jacobi polynomial, i.e.

$ P_n^{\alpha,\beta}(z) = 0 $

This routine is only value for $( \alpha > -1, \beta > -1)$

and uses polynomial deflation in a Newton iteration

Here is the call graph for this function:

LIB_UTILITIES_EXPORT void Polylib::JacZeros ( const int  n,
double *  a,
double *  b,
const double  alpha,
const double  beta 
)

Zero and Weight determination through the eigenvalues and eigenvectors of a tridiagonal.

matrix from the three term recurrence relationship.

Set up a symmetric tridiagonal matrix

$ \left [ \begin{array}{ccccc} a[0] & b[0] & & & \\ b[0] & a[1] & b[1] & & \\ 0 & \ddots & \ddots & \ddots & \\ & & \ddots & \ddots & b[n-2] \\ & & & b[n-2] & a[n-1] \end{array} \right ] $

Where the coefficients a[n], b[n] come from the recurrence relation

$ b_j p_j(z) = (z - a_j ) p_{j-1}(z) - b_{j-1} p_{j-2}(z) $

where $ j=n+1$ and $p_j(z)$ are the Jacobi (normalized)

orthogonal polynomials $ \alpha,\beta > -1$( integer values and

halves). Since the polynomials are orthonormalized, the tridiagonal

matrix is guaranteed to be symmetric. The eigenvalues of this

matrix are the zeros of the Jacobi polynomial.

Here is the call graph for this function:

void Polylib::JKMatrix ( int  n,
double *  a,
double *  b 
)

Calcualtes the Jacobi-kronrod matrix by determining the.

a and coefficients.

The first \a 3n+1 coefficients are already known



For more information refer to:

"Dirk P. Laurie, Calcualtion of Gauss-Kronrod quadrature rules"

Here is the caller graph for this function:

static void Polylib::RecCoeff ( const int  n,
double *  a,
double *  b,
const double  alpha,
const double  beta 
)
static

The routine finds the recurrence coefficients a and.

b of the orthogonal polynomials

Here is the call graph for this function:

Here is the caller graph for this function:

static void Polylib::TriQL ( const int  n,
double *  d,
double *  e,
double **  z 
)
static

QL algorithm for symmetric tridiagonal matrix.

This subroutine is a translation of an algol procedure,

num. math. 12, 377-383(1968) by martin and wilkinson, as modified

in num. math. 15, 450(1970) by dubrulle. Handbook for

auto. comp., vol.ii-linear algebra, 241-248(1971). This is a

modified version from numerical recipes.

This subroutine finds the eigenvalues and first components of the

eigenvectors of a symmetric tridiagonal matrix by the implicit QL

method.

on input:

  • n is the order of the matrix;
  • d contains the diagonal elements of the input matrix;
  • e contains the subdiagonal elements of the input matrix

in its first n-2 positions.

- z is the n by n identity matrix

on output:

  • d contains the eigenvalues in ascending order.
  • e contains the weight values - modifications of the first component
    of normalised eigenvectors

Here is the caller graph for this function:

LIB_UTILITIES_EXPORT void Polylib::zwgj ( double *  z,
double *  w,
const int  np,
const double  alpha,
const double  beta 
)

Gauss-Jacobi zeros and weights.

  • Generate np Gauss Jacobi zeros, z, and weights,w,

associated with the Jacobi polynomial $ P^{\alpha,\beta}_{np}(z) $,

  • Exact for polynomials of order 2np-1 or less

Here is the call graph for this function:

LIB_UTILITIES_EXPORT void Polylib::zwgk ( double *  z,
double *  w,
const int  npt,
const double  alpha,
const double  beta 
)

Gauss-Kronrod-Jacobi zeros and weights.

  • Generate npt=2*np+1 Gauss-Kronrod Jacobi zeros, z, and weights,w,

associated with the Jacobi polynomial $ P^{\alpha,\beta}_{np}(z) $,

  • Exact for polynomials of order 3np+1 or less

Here is the call graph for this function:

LIB_UTILITIES_EXPORT void Polylib::zwglj ( double *  z,
double *  w,
const int  np,
const double  alpha,
const double  beta 
)

Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1.

  • Generate np Gauss-Lobatto-Jacobi points, z, and weights, w,

associated with polynomial $ (1-z)(1+z) P^{\alpha+1,\beta+1}_{np-2}(z) $

  • Exact for polynomials of order 2np-3 or less

Here is the call graph for this function:

LIB_UTILITIES_EXPORT void Polylib::zwgrjm ( double *  z,
double *  w,
const int  np,
const double  alpha,
const double  beta 
)

Gauss-Radau-Jacobi zeros and weights with end point at z=-1.

  • Generate np Gauss-Radau-Jacobi zeros, z, and weights,w,

associated with the polynomial $(1+z) P^{\alpha,\beta+1}_{np-1}(z) $.

  • Exact for polynomials of order 2np-2 or less

Here is the call graph for this function:

LIB_UTILITIES_EXPORT void Polylib::zwgrjp ( double *  z,
double *  w,
const int  np,
const double  alpha,
const double  beta 
)

Gauss-Radau-Jacobi zeros and weights with end point at z=1.

  • Generate np Gauss-Radau-Jacobi zeros, z, and weights,w,

associated with the polynomial $(1-z) P^{\alpha+1,\beta}_{np-1}(z) $.

  • Exact for polynomials of order 2np-2 or less

Here is the call graph for this function:

LIB_UTILITIES_EXPORT void Polylib::zwlk ( double *  z,
double *  w,
const int  npt,
const double  alpha,
const double  beta 
)

Gauss-Lobatto-Kronrod-Jacobi zeros and weights.

  • Generate npt=2*np-1 Lobatto-Kronrod Jacobi zeros, z, and weights,w,

associated with the Jacobi polynomial $ P^{\alpha,\beta}_{np}(z) $,

Here is the call graph for this function:

LIB_UTILITIES_EXPORT void Polylib::zwrk ( double *  z,
double *  w,
const int  npt,
const double  alpha,
const double  beta 
)

Gauss-Radau-Kronrod-Jacobi zeros and weights.

  • Generate npt=2*np Radau-Kronrod Jacobi zeros, z, and weights,w,

associated with the Jacobi polynomial $ P^{\alpha,\beta}_{np}(z) $,

Here is the call graph for this function: