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.
static void	TriQL (const int n, double d, double e, double **z)
	QL algorithm for symmetric tridiagonal matrix.
double	gammaF (const double x)
	Calculate the Gamma function , $\Gamma(n)$ , for integer.
static void	RecCoeff (const int n, double a, double b, const double alpha, const double beta)
	The routine finds the recurrence coefficients a and.
void	JKMatrix (int n, double a, double b)
	Calcualtes the Jacobi-kronrod matrix by determining the.
void	chri1 (int n, double a, double b, double a0, double b0, double z)
	Given a weight function through the first n+1.
void	zwgj (double z, double w, const int np, const double alpha, const double beta)
	Gauss-Jacobi zeros and weights.
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.
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.
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.
void	zwgk (double z, double w, const int npt, const double alpha, const double beta)
	Gauss-Kronrod-Jacobi zeros and weights.
void	zwrk (double z, double w, const int npt, const double alpha, const double beta)
	Gauss-Radau-Kronrod-Jacobi zeros and weights.
void	zwlk (double z, double w, const int npt, const double alpha, const double beta)
	Gauss-Lobatto-Kronrod-Jacobi zeros and weights.
void	Dgj (double D, const double z, const int np, const double alpha, const double beta)
	Compute the Derivative Matrix and its transpose associated.
void	Dgrjm (double D, const double z, const int np, const double alpha, const double beta)
	Compute the Derivative Matrix and its transpose associated.
void	Dgrjp (double D, const double z, const int np, const double alpha, const double beta)
	Compute the Derivative Matrix associated with the.
void	Dglj (double D, const double z, const int np, const double alpha, const double beta)
	Compute the Derivative Matrix associated with the.
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.
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.
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.
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.
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.
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.
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.
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.
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)$ .
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.
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.

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.

Definition at line 2926 of file Polylib.cpp.

Referenced by zwlk(), and zwrk().

    {
        
        double q = ceil(3.0*n/2);
        int size = (int)q+1;
        if(size < n+1)
        {
            fprintf(stderr,"input arrays a and b are too short\n");
        }
        double* r = new double[n+1];
        r[0] = z - a[0];
        r[1] = z - a[1] - b[1]/r[0];
        a0[0] = a[1] + r[1] - r[0];
        b0[0] = -r[0]*b[0];
        if(n == 1)
        {
            return;
        }
        int k = 0;
        for(k = 1; k < n; k++)
        {
            r[k+1] = z - a[k+1] - b[k+1]/r[k];
            a0[k] = a[k+1] + r[k+1] - r[k];
            b0[k] = b[k] * r[k]/r[k-1];
        }
        
    }

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])$

Definition at line 932 of file Polylib.cpp.

References jacobd().

Referenced by Nektar::LibUtilities::GaussPoints::CalculateDerivMatrix(), and main().

    {
        double one = 1.0, two = 2.0;
        if (np <= 0){
            D[0] = 0.0;
        }
        else{
            register int i,j; 
            double *pd;
            pd = (double *)malloc(np*sizeof(double));
            jacobd(np,z,pd,np,alpha,beta);
            for (i = 0; i < np; i++){
                for (j = 0; j < np; j++){
                    if (i != j) 
                        D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
                    else    
                        D[i*np+j] = (alpha - beta + (alpha + beta + two)*z[j])/
                        (two*(one - z[j]*z[j]));
                }
            }
            free(pd);
        }
        return;
    }

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])$

Definition at line 1216 of file Polylib.cpp.

References gammaF(), and jacobd().

Referenced by Nektar::LibUtilities::GaussPoints::CalculateDerivMatrix(), and main().

    {
        if (np <= 0){
            D[0] = 0.0;
        }
        else{
            register int i, j; 
            double   one = 1.0, two = 2.0;
            double   *pd;
            pd  = (double *)malloc(np*sizeof(double));
            pd[0]  = two*pow(-one,np)*gammaF(np + beta);
            pd[0] /= gammaF(np - one)*gammaF(beta + two);
            jacobd(np-2,z+1,pd+1,np-2,alpha+1,beta+1);
            for(i = 1; i < np-1; ++i) pd[i] *= (one-z[i]*z[i]);
            pd[np-1]  = -two*gammaF(np + alpha);
            pd[np-1] /= gammaF(np - one)*gammaF(alpha + two);
            for (i = 0; i < np; i++)
                for (j = 0; j < np; j++){
                    if (i != j) 
                        D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
                    else { 
                        if (j == 0)
                            D[i*np+j] = (alpha - (np-1)*(np + alpha + beta))/(two*(beta+ two));
                        else if (j == np-1)
                            D[i*np+j] =-(beta - (np-1)*(np + alpha + beta))/(two*(alpha+ two));
                        else
                            D[i*np+j] = (alpha - beta + (alpha + beta)*z[j])/
                            (two*(one - z[j]*z[j]));
                    }
                }
                free(pd);
        }
        return;
    }

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])$

Definition at line 1018 of file Polylib.cpp.

References gammaF(), and jacobd().

Referenced by Nektar::LibUtilities::GaussPoints::CalculateDerivMatrix(), and main().

    {
        if (np <= 0){
            D[0] = 0.0;
        }
        else{
            register int i, j; 
            double   one = 1.0, two = 2.0;
            double   *pd;
            pd  = (double *)malloc(np*sizeof(double));
            pd[0] = pow(-one,np-1)*gammaF(np+beta+one);
            pd[0] /= gammaF(np)*gammaF(beta+two);
            jacobd(np-1,z+1,pd+1,np-1,alpha,beta+1);
            for(i = 1; i < np; ++i) pd[i] *= (1+z[i]);
            for (i = 0; i < np; i++)
                for (j = 0; j < np; j++){
                    if (i != j) 
                        D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
                    else { 
                        if(j == 0)
                            D[i*np+j] = -(np + alpha + beta + one)*(np - one)/
                            (two*(beta + two));
                        else
                            D[i*np+j] = (alpha - beta + one + (alpha + beta + one)*z[j])/
                            (two*(one - z[j]*z[j]));
                    }
                }
                free(pd);
        }
        return;
    }

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])$

Definition at line 1116 of file Polylib.cpp.

References gammaF(), and jacobd().

Referenced by Nektar::LibUtilities::GaussPoints::CalculateDerivMatrix(), and main().

    {
        if (np <= 0){
            D[0] = 0.0;
        }
        else{
            register int i, j; 
            double   one = 1.0, two = 2.0;
            double   *pd;
            pd  = (double *)malloc(np*sizeof(double));
            jacobd(np-1,z,pd,np-1,alpha+1,beta);
            for(i = 0; i < np-1; ++i) pd[i] *= (1-z[i]);
            pd[np-1] = -gammaF(np+alpha+one);
            pd[np-1] /= gammaF(np)*gammaF(alpha+two);
            for (i = 0; i < np; i++)
                for (j = 0; j < np; j++){
                    if (i != j) 
                        D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
                    else { 
                        if(j == np-1)
                            D[i*np+j] = (np + alpha + beta + one)*(np - one)/
                            (two*(alpha + two));
                        else
                            D[i*np+j] = (alpha - beta - one + (alpha + beta + one)*z[j])/
                            (two*(one - z[j]*z[j]));
                    }
                }
                free(pd);
        }
        return;
    }

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)$

Definition at line 2200 of file Polylib.cpp.

Referenced by Dglj(), Dgrjm(), Dgrjp(), RecCoeff(), zwgj(), zwglj(), zwgrjm(), and zwgrjp().

                                 {
        double gamma = 1.0;
        if     (x == -0.5) gamma = -2.0*sqrt(M_PI);
        else if (!x) return gamma;
        else if ((x-(int)x) == 0.5){ 
            int n = (int) x;
            double tmp = x;
            gamma = sqrt(M_PI);
            while(n--){
                tmp   -= 1.0;
                gamma *= tmp;
            }
        }
        else if ((x-(int)x) == 0.0){
            int n = (int) x;
            double tmp = x;
            while(--n){
                tmp   -= 1.0;
                gamma *= tmp;
            }
        }  
        else
            fprintf(stderr,"%lf is not of integer or half order\n",x);
        return gamma;
    }

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}$

Definition at line 1340 of file Polylib.cpp.

References EPS, jacobd(), and jacobfd().

Referenced by Nektar::LibUtilities::Basis::GenBasis(), and Imgj().

    {
        double zi, dz, p, pd, h;
        zi  = *(zgj+i);
        dz  = z - zi;
        if (fabs(dz) < EPS) return 1.0;
        jacobd (1, &zi, &pd , np, alpha, beta);
        jacobfd(1, &z , &p, NULL , np, alpha, beta);
        h = p/(pd*dz);
        return h;
    }

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}$

Definition at line 1584 of file Polylib.cpp.

References EPS, jacobd(), and jacobfd().

Referenced by Nektar::LibUtilities::Basis::GenBasis(), and Imglj().

    {
        double one = 1., two = 2.;
        double zi, dz, p, pd, h;
        zi  = *(zglj+i);
        dz  = z - zi;
        if (fabs(dz) < EPS) return 1.0;
        jacobfd(1, &zi, &p , NULL, np-2, alpha + one, beta + one);
        // need to use this routine in caes z = -1 or 1
        jacobd (1, &zi, &pd, np-2, alpha + one, beta + one);
        h = (one - zi*zi)*pd - two*zi*p;
        jacobfd(1, &z, &p, NULL, np-2, alpha + one, beta + one);
        h = (one - z*z)*p/(h*dz);
        return h;
    }

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}$

Definition at line 1416 of file Polylib.cpp.

References EPS, jacobd(), and jacobfd().

Referenced by Imgrjm().

    {
        double zi, dz, p, pd, h;
        zi  = *(zgrj+i);
        dz  = z - zi;
        if (fabs(dz) < EPS) return 1.0;
        jacobfd (1, &zi, &p , NULL, np-1, alpha, beta + 1);
        // need to use this routine in caes zi = -1 or 1
        jacobd  (1, &zi, &pd, np-1, alpha, beta + 1);
        h = (1.0 + zi)*pd + p;
        jacobfd (1, &z, &p, NULL,  np-1, alpha, beta + 1);
        h = (1.0 + z )*p/(h*dz);
        return h;
    }

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}$

Definition at line 1500 of file Polylib.cpp.

References EPS, jacobd(), and jacobfd().

Referenced by Imgrjp().

    {
        double zi, dz, p, pd, h;
        zi  = *(zgrj+i);
        dz  = z - zi;
        if (fabs(dz) < EPS) return 1.0;
        jacobfd (1, &zi, &p , NULL, np-1, alpha+1, beta );
        // need to use this routine in caes z = -1 or 1
        jacobd  (1, &zi, &pd, np-1, alpha+1, beta );
        h = (1.0 - zi)*pd - p;
        jacobfd (1, &z, &p, NULL,  np-1, alpha+1, beta);
        h = (1.0 - z )*p/(h*dz);
        return h;
    }

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])$

Definition at line 1652 of file Polylib.cpp.

References hgj().

Referenced by Nektar::LibUtilities::GaussPoints::CalculateInterpMatrix(), and main().

                                                           {
            double zp;
            register int i, j;
            for (i = 0; i < nz; ++i) {
                for (j = 0; j < mz; ++j)
                {
                    zp = zm[j];
                    im [i*mz+j] = hgj(i, zp, zgj, nz, alpha, beta);
                }
            }
            return;
    }

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])$

Definition at line 1832 of file Polylib.cpp.

References hglj().

Referenced by Nektar::LibUtilities::GaussPoints::CalculateInterpMatrix(), and main().

    {
        double zp;
        register int i, j;
        for (i = 0; i < nz; i++) {
            for (j = 0; j < mz; j++)
            {
                zp = zm[j];
                im[i*mz+j] = hglj(i, zp, zglj, nz, alpha, beta);
            }
        }
        return;
    }

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])$

Definition at line 1710 of file Polylib.cpp.

References hgrjm().

Referenced by Nektar::LibUtilities::GaussPoints::CalculateInterpMatrix(), and main().

    {
        double zp;
        register int i, j;
        for (i = 0; i < nz; i++) {
            for (j = 0; j < mz; j++)
            {
                zp = zm[j];
                im [i*mz+j] = hgrjm(i, zp, zgrj, nz, alpha, beta);
            }
        }
        return;
    }

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])$

Definition at line 1770 of file Polylib.cpp.

References hgrjp().

Referenced by Nektar::LibUtilities::GaussPoints::CalculateInterpMatrix(), and main().

    {
            double zp;
            register int i, j;
            for (i = 0; i < nz; i++) {
                for (j = 0; j < mz; j++)
                {
                    zp = zm[j];
                    im [i*mz+j] = hgrjp(i, zp, zgrj, nz, alpha, beta);
                }
            }
            return;
    }

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$

Definition at line 2146 of file Polylib.cpp.

References jacobfd().

Referenced by Dgj(), Dglj(), Dgrjm(), Dgrjp(), hgj(), hglj(), hgrjm(), hgrjp(), Nektar::SolverUtils::AdvectionFR::v_SetupCFunctions(), Nektar::SolverUtils::DiffusionLFR::v_SetupCFunctions(), Nektar::SolverUtils::DiffusionLFRNS::v_SetupCFunctions(), and zwgj().

    {
        register int i;
        double one = 1.0;
        if(n == 0)
            for(i = 0; i < np; ++i) polyd[i] = 0.0;
        else{
            //jacobf(np,z,polyd,n-1,alpha+one,beta+one);
            jacobfd(np,z,polyd,NULL,n-1,alpha+one,beta+one);
            for(i = 0; i < np; ++i) polyd[i] *= 0.5*(alpha + beta + (double)n + one);
        }
        return;
    }

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.

Definition at line 1946 of file Polylib.cpp.

Referenced by Nektar::LibUtilities::Basis::GenBasis(), hgj(), hglj(), hgrjm(), hgrjp(), jacobd(), Jacobz(), main(), Nektar::MultiRegions::ExpList1D::PeriodicEval(), zwglj(), zwgrjm(), and zwgrjp().

                                                           {
            register int i;
            double  zero = 0.0, one = 1.0, two = 2.0;
            if(!np)
                return;
            if(n == 0){
                if(poly_in)
                    for(i = 0; i < np; ++i) 
                        poly_in[i] = one;
                if(polyd)
                    for(i = 0; i < np; ++i) 
                        polyd[i] = zero; 
            }
            else if (n == 1){
                if(poly_in)
                    for(i = 0; i < np; ++i) 
                        poly_in[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
                if(polyd)
                    for(i = 0; i < np; ++i) 
                        polyd[i] = 0.5*(alpha + beta + two);
            }
            else{
                register int k;
                double   a1,a2,a3,a4;
                double   two = 2.0, apb = alpha + beta;
                double   *poly, *polyn1,*polyn2;
                if(poly_in){ // switch for case of no poynomial function return
                    polyn1 = (double *)malloc(2*np*sizeof(double));
                    polyn2 = polyn1+np; 
                    poly   = poly_in;
                }
                else{
                    polyn1 = (double *)malloc(3*np*sizeof(double));
                    polyn2 = polyn1+np; 
                    poly   = polyn2+np;      
                }
                for(i = 0; i < np; ++i){
                    polyn2[i] = one;
                    polyn1[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
                }
                for(k = 2; k <= n; ++k){
                    a1 =  two*k*(k + apb)*(two*k + apb - two);
                    a2 = (two*k + apb - one)*(alpha*alpha - beta*beta);
                    a3 = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb);
                    a4 =  two*(k + alpha - one)*(k + beta - one)*(two*k + apb);
                    a2 /= a1;
                    a3 /= a1;
                    a4 /= a1;
                    for(i = 0; i < np; ++i){
                        poly  [i] = (a2 + a3*z[i])*polyn1[i] - a4*polyn2[i];
                        polyn2[i] = polyn1[i];
                        polyn1[i] = poly  [i];
                    }
                }
                if(polyd){
                    a1 = n*(alpha - beta);
                    a2 = n*(two*n + alpha + beta);
                    a3 = two*(n + alpha)*(n + beta);
                    a4 = (two*n + alpha + beta);
                    a1 /= a4;  a2 /= a4;   a3 /= a4;
                    // note polyn2 points to polyn1 at end of poly iterations
                    for(i = 0; i < np; ++i){
                        polyd[i]  = (a1- a2*z[i])*poly[i] + a3*polyn2[i];
                        polyd[i] /= (one - z[i]*z[i]);
                    }
                }
                free(polyn1);
            }
            return;
    }

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

Definition at line 2274 of file Polylib.cpp.

References EPS, jacobfd(), and STOP.

                          {
            register int i,j,k;
            double   dth = M_PI/(2.0*(double)n);
            double   poly,pder,rlast=0.0;
            double   sum,delr,r;
            double one = 1.0, two = 2.0;
            if(!n)
                return;
            for(k = 0; k < n; ++k){
                r = -cos((two*(double)k + one) * dth);
                if(k) r = 0.5*(r + rlast);
                for(j = 1; j < STOP; ++j){
                    jacobfd(1,&r,&poly, &pder, n, alpha, beta);
                    for(i = 0, sum = 0.0; i < k; ++i) sum += one/(r - z[i]);
                    delr = -poly / (pder - sum * poly);
                    r   += delr;
                    if( fabs(delr) < EPS ) break;
                }
                z[k]  = r;
                rlast = r;
            }
            return;
    }

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.

Definition at line 2384 of file Polylib.cpp.

References RecCoeff(), and TriQL().

                          {
            
            int i,j;
            RecCoeff(n,a,b,alpha,beta);
            
            double **z = new double*[n];
            for(i = 0; i < n; i++)
            {
                z[i] = new double[n];
                for(j = 0; j < n; j++)
                {
                    z[i][j] = 0.0;
                }
            }
            for(i = 0; i < n; i++)
            {
                z[i][i] = 1.0;
            }
                    
            // find eigenvalues and eigenvectors
            TriQL(n, a, b,z);
            return;
    }

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

Calcualtes the Jacobi-kronrod matrix by determining the.

a and coefficients.

The first 3n+1 coefficients are already known

For more information refer to:

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

Definition at line 2764 of file Polylib.cpp.

Referenced by zwgk(), zwlk(), and zwrk().

    {
        int i,j,k,m;
        // Working storage
        int size = (int)floor(n/2.0)+2;
        double *s = new double[size];
        double *t = new double[size];
        // Initialize s and t to zero
        for(i = 0; i < size; i++)
        {
            s[i] = 0.0;
            t[i] = 0.0;
        }
        t[1] = b[n+1];
        for(m = 0; m <= n-2; m++)
        {
             double u = 0.0;
             for(k = (int)floor((m+1)/2.0); k >= 0; k--)
             {
                int l = m-k;
                u = u+(a[k+n+1]-a[l])*t[k+1] + b[k+n+1]*s[k] - b[l]*s[k+1];
                s[k+1] = u;
             }
             // Swap the contents of s and t
             double *hold = s;
             s = t;
             t  = hold;
        }
        
        
        for(j = (int)floor(n/2.0); j >= 0; j--)
        {
            s[j+1] = s[j];
        }
        for(m = n-1; m <= 2*n-3; m++)
        {
            double u = 0;
            for(k = m+1-n; k <= floor((m-1)/2.0); k++)
            {
                int l = m-k;
                j = n-1-l;
                u = u-(a[k+n+1]-a[l])*t[j+1] - b[k+n+1]*s[j+1] + b[l]*s[j+2];
                s[j+1] = u;
            }
            if(m%2 == 0)
            {
                k = m/2;
                a[k+n+1] = a[k] + (s[j+1]-b[k+n+1]*s[j+2])/t[j+2];
                
            }else
            {
                k = (m+1)/2;
                b[k+n+1] = s[j+1]/s[j+2];
            }
            
            
            // Swap the contents of s and t
            double  *hold = s;
            s = t;
            t  = hold;
        }
        a[2*n ] = a[n-1]-b[2*n]*s[1]/t[1];
        
    }

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

Definition at line 2444 of file Polylib.cpp.

References gammaF().

Referenced by JacZeros(), zwgk(), zwlk(), and zwrk().

                                {
        int i;
        double apb, apbi,a2b2;
        
        if(!n)
            return;
        // generate normalised terms 
        apb  = alpha + beta;
        apbi = 2.0 + apb;
        b[0] = pow(2.0,apb+1.0)*gammaF(alpha+1.0)*gammaF(beta+1.0)/gammaF(apbi); //MuZero
        a[0]   = (beta-alpha)/apbi;
        b[1]   = (4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));
        
        a2b2 = beta*beta-alpha*alpha;
        
        for(i = 1; i < n-1; i++){
            apbi = 2.0*(i+1) + apb;
            a[i] = a2b2/((apbi-2.0)*apbi);
            b[i+1] = (4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
                 ((apbi*apbi-1)*apbi*apbi));
        }
        
        apbi   = 2.0*n + apb;
        a[n-1] = a2b2/((apbi-2.0)*apbi);
        
    }

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

Definition at line 2560 of file Polylib.cpp.

References sign, and STOP.

Referenced by JacZeros(), zwgk(), zwlk(), and zwrk().

                                                                   {
        int m,l,iter,i,k;
        double s,r,p,g,f,dd,c,b;
        
        double MuZero = e[0];
        
        // Renumber the elements of e
        for(i = 0; i < n-1; i++)
        {
            e[i] = sqrt(e[i+1]);
        }
        e[n-1] = 0.0;
        for (l=0;l<n;l++) {
            iter=0;
            do {
                for (m=l;m<n-1;m++) {
                    dd=fabs(d[m])+fabs(d[m+1]);
                    if (fabs(e[m])+dd == dd) break;
                }
                if (m != l) {
                    if (iter++ == STOP){
                        fprintf(stderr,"triQL: Too many iterations in TQLI");
                        exit(1);
                    }
                    g=(d[l+1]-d[l])/(2.0*e[l]);
                    r=sqrt((g*g)+1.0);
                    g=d[m]-d[l]+e[l]/(g+sign(r,g));
                    s=c=1.0;
                    p=0.0;
                    for (i=m-1;i>=l;i--) {
                        f=s*e[i];
                        b=c*e[i];
                        if (fabs(f) >= fabs(g)) {
                            c=g/f;
                            r=sqrt((c*c)+1.0);
                            e[i+1]=f*r;
                            c *= (s=1.0/r);
                        } else {
                            s=f/g;
                            r=sqrt((s*s)+1.0);
                            e[i+1]=g*r;
                            s *= (c=1.0/r);
                        }
                        g=d[i+1]-p;
                        r=(d[i]-g)*s+2.0*c*b;
                        p=s*r;
                        d[i+1]=g+p;
                        g=c*r-b;
                        // Calculate the eigenvectors
                        for(k = 0; k < n; k++)
                        {
                            f = z[k][i+1];
                            z[k][i+1] = s*z[k][i] + c*f;
                            z[k][i] = c*z[k][i] - s*f;
                        }
                    
                    }
                    d[l]=d[l]-p;
                    e[l]=g;
                    e[m]=0.0;
                }
            } while (m != l);
        }
        // order eigenvalues
        // Since we only need the first component of the eigenvectors
        // to calcualte the weight, we only swap the first components
        for(i = 0; i < n-1; ++i){ 
            k = i;
            p = d[i];
            for(l = i+1; l < n; ++l)
                if (d[l] < p) {
                    k = l;
                    p = d[l];
                }
            d[k] = d[i]; 
            d[i] = p;
            double temp = z[0][k];
            z[0][k] = z[0][i];
            z[0][i] = temp;
        }
        // Calculate the weights
        for(i =0 ; i < n; i++)
        {
            e[i] = MuZero*z[0][i]*z[0][i];
        }
    }

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

Definition at line 102 of file Polylib.cpp.

References gammaF(), jacobd(), and jacobz.

Referenced by Nektar::LibUtilities::GaussPoints::CalculatePoints(), and main().

    {
        register int i;
        double fac, one = 1.0, two = 2.0, apb = alpha + beta;
        jacobz (np,z,alpha,beta);
        jacobd (np,z,w,np,alpha,beta);
        fac  = pow(two,apb + one)*gammaF(alpha + np + one)*gammaF(beta + np + one);
        fac /= gammaF(np + one)*gammaF(apb + np + one);
        for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]*(one-z[i]*z[i]));
        return;
    }

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

Definition at line 380 of file Polylib.cpp.

References JKMatrix(), RecCoeff(), and TriQL().

Referenced by Nektar::LibUtilities::GaussPoints::CalculatePoints().

    {
        
        int np = (npt-1)/2; 
        
        int  i,j;
        
        // number of kronrod points associated with the np gauss rule
        int kpoints = 2*np + 1;
        // Define the number of required recurrence coefficents
        int ncoeffs = (int)floor(3.0*(np+1)/2);
        // Define arrays  for the recurrence coefficients
        // We will use these arrays for the Kronrod results too, hence the 
        // reason for the size of the arrays
        double *a = new double[kpoints];
        double *b = new double[kpoints];
        // Initialize a and b to zero
        for(i = 0; i < kpoints; i++)
        {
            a[i] = 0.0;
            b[i] = 0.0;
        }
        // Call the routine to calculate the recurrence coefficients
        RecCoeff(ncoeffs,a,b,alpha,beta);
        // Call the routine to caluclate the jacobi-Kronrod matrix
        JKMatrix(np,a,b);
        // Set up the identity matrix
        double** zmatrix = new double*[kpoints];
        for(i = 0; i < kpoints; i++)
        {
            zmatrix[i] = new double[kpoints];
            for(j = 0; j < kpoints; j++)
            {
                zmatrix[i][j] = 0.0;
            }
        }
        for(i = 0; i < kpoints; i++)
        {
            zmatrix[i][i] = 1.0;
        }
        // Calculte the points and weights
        TriQL(kpoints, a, b, zmatrix);
        for(i = 0; i < kpoints; i++)
        {
            z[i] = a[i];
            w[i] = b[i];
        }
        
    }

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

Definition at line 306 of file Polylib.cpp.

References gammaF(), jacobfd(), and jacobz.

Referenced by Nektar::LibUtilities::GaussPoints::CalculatePoints(), and main().

    {
        if( np == 1 ){
            z[0] = 0.0;
            w[0] = 2.0;
        }
        else{
            register int i;
            double   fac, one = 1.0, apb = alpha + beta, two = 2.0;
            z[0]    = -one;
            z[np-1] =  one;
            jacobz  (np-2,z + 1,alpha + one,beta + one); 
            jacobfd (np,z,w,NULL,np-1,alpha,beta);
            fac  = pow(two,apb + 1)*gammaF(alpha + np)*gammaF(beta + np);
            fac /= (np-1)*gammaF(np)*gammaF(alpha + beta + np + one);
            for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]);
            w[0]    *= (beta  + one);
            w[np-1] *= (alpha + one);
        }
        return;
    }

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

Definition at line 158 of file Polylib.cpp.

References gammaF(), jacobfd(), and jacobz.

Referenced by Nektar::LibUtilities::GaussPoints::CalculatePoints(), and main().

    {
        if(np == 1){
            z[0] = 0.0;
            w[0] = 2.0;
        }
        else{
            register int i;
            double fac, one = 1.0, two = 2.0, apb = alpha + beta;
            z[0] = -one;
            jacobz  (np-1,z+1,alpha,beta+1);
            jacobfd (np,z,w,NULL,np-1,alpha,beta);
            fac  = pow(two,apb)*gammaF(alpha + np)*gammaF(beta + np);
            fac /= gammaF(np)*(beta + np)*gammaF(apb + np + 1);
            for(i = 0; i < np; ++i) w[i] = fac*(1-z[i])/(w[i]*w[i]);
            w[0] *= (beta + one);
        }
        return;
    }

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

Definition at line 234 of file Polylib.cpp.

References gammaF(), jacobfd(), and jacobz.

Referenced by Nektar::LibUtilities::GaussPoints::CalculatePoints(), and main().

    {
        if(np == 1){
            z[0] = 0.0;
            w[0] = 2.0;
        }
        else{
            register int i;
            double fac, one = 1.0, two = 2.0, apb = alpha + beta;
            jacobz  (np-1,z,alpha+1,beta);
            z[np-1] = one;
            jacobfd (np,z,w,NULL,np-1,alpha,beta);
            fac  = pow(two,apb)*gammaF(alpha + np)*gammaF(beta + np);
            fac /= gammaF(np)*(alpha + np)*gammaF(apb + np + 1);
            for(i = 0; i < np; ++i) w[i] = fac*(1+z[i])/(w[i]*w[i]);
            w[np-1] *= (alpha + one);
        }
        return;
    }

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)$ ,

Definition at line 706 of file Polylib.cpp.

References chri1(), JKMatrix(), RecCoeff(), and TriQL().

Referenced by Nektar::LibUtilities::GaussPoints::CalculatePoints().

    {
        int np = (npt+1)/2;
        
        if(np < 4)
        {
            fprintf (stderr,"too few points in formula\n");
            return;
        }
        
        double endl = -1;
        double endr = 1;
        int i,j;
        // number of kronrod points associated with the np gauss rule
        int kpoints = 2*np-1;
        // Define the number of required recurrence coefficents
        int ncoeffs = (int)ceil(3.0*np/2)-1;
        // Define arrays  for the recurrence coefficients
        double *a = new double[ncoeffs+1];
        double *b = new double[ncoeffs+1];
        // Initialize a and b to zero
        for(i = 0; i < ncoeffs+1; i++)
        {
            a[i] = 0.0;
            b[i] = 0.0;
        }
        // Call the routine to calculate the recurrence coefficients
        RecCoeff(ncoeffs,a,b,alpha,beta);
        double* a0 = new double[ncoeffs];
        double* b0 = new double[ncoeffs];
        chri1(ncoeffs,a,b,a0,b0,endl);
        double* a1 = new double[ncoeffs-1];
        double* b1 = new double[ncoeffs-1];
        chri1(ncoeffs-1,a0,b0,a1,b1,endr);
        double s = b1[0]/fabs(b1[0]);
        b1[0] = s*b1[0];
        // Finding the 2*np-1 gauss-kronrod points
        double* z1 = new double[2*np-3];
        double* w1 = new double[2*np-3];
        for(i = 0; i < ncoeffs; i++)
        {
            z1[i] = a1[i];
            w1[i] = b1[i];
        }
        JKMatrix(np-2,z1,w1);
        // Set up the identity matrix
        double** zmatrix = new double*[2*np-3];
        for(i = 0; i < 2*np-3; i++)
        {
            zmatrix[i] = new double[2*np-3];
            for(j = 0; j < 2*np-3; j++)
            {
                zmatrix[i][j] = 0.0;
            }
        }
        for(i = 0; i < 2*np-3; i++)
        {
            zmatrix[i][i] = 1.0;
        }
        // Calculate the points and weights
        TriQL(2*np-3, z1, w1, zmatrix);
        double sumW1 = 0.0;
        double sumW1Z1 = 0.0;
        for(i = 0; i < 2*np-3; i++)
        {
            w1[i] = s*w1[i]/(z1[i]-endl)/(z1[i]-endr);
            sumW1 += w1[i];
            sumW1Z1 += z1[i]*w1[i];
        }
        double c0 = b[0]-sumW1;
        double c1 = a[0]*b[0]-sumW1Z1;
            
        z[0] = endl;
        z[2*np-2] = endr;
        w[0] = (c0*endr-c1)/(endr-endl);
        w[2*np-2] = (c1-c0*endl)/(endr-endl); 
        for(i = 1; i < kpoints-1; i++)
        {
            z[i] = z1[i-1];
            w[i] = w1[i-1];
        }
    }

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)$ ,

Definition at line 514 of file Polylib.cpp.

References chri1(), JKMatrix(), RecCoeff(), and TriQL().

Referenced by Nektar::LibUtilities::GaussPoints::CalculatePoints().

    {
        
        int np = npt/2;
        if(np < 2)
        {
            fprintf(stderr,"too few points in formula\n");
            return;
        }
        double end0 = -1;
        int i,j;
        // number of kronrod points associated with the np gauss rule
        int kpoints = 2*np;
        // Define the number of required recurrence coefficents
        int ncoeffs = (int)ceil(3.0*np/2);
        // Define arrays  for the recurrence coefficients
        double *a = new double[ncoeffs+1];
        double *b = new double[ncoeffs+1];
        // Initialize a and b to zero
        for(i = 0; i < ncoeffs+1; i++)
        {
            a[i] = 0.0;
            b[i] = 0.0;
        }
        // Call the routine to calculate the recurrence coefficients
        RecCoeff(ncoeffs,a,b,alpha,beta);
        double* a0 = new double[ncoeffs];
        double* b0 = new double[ncoeffs];
        chri1(ncoeffs,a,b,a0,b0,end0);
        double s = b0[0]/fabs(b0[0]);
        b0[0] = s*b0[0];
        // Finding the 2*np-1 gauss-kronrod points
        double* z1 = new double[2*np-1];
        double* w1 = new double[2*np-1];
        for(i = 0; i < ncoeffs; i++)
        {
            z1[i] = a0[i];
            w1[i] = b0[i];
        }
        JKMatrix(np-1,z1,w1);
        // Set up the identity matrix
        double** zmatrix = new double*[2*np-1];
        for(i = 0; i < 2*np-1; i++)
        {
            zmatrix[i] = new double[2*np-1];
            for(j = 0; j < 2*np-1; j++)
            {
                zmatrix[i][j] = 0.0;
            }
        }
        for(i = 0; i < 2*np-1; i++)
        {
            zmatrix[i][i] = 1.0;
        }
        // Calculate the points and weights
        TriQL(2*np-1, z1, w1, zmatrix);
        double sumW1 = 0.0;
        for(i = 0; i < 2*np-1; i++)
        {
            w1[i] = s*w1[i]/(z1[i]-end0);
            sumW1 += w1[i];
        }
            
        z[0] = end0;
        w[0] = b[0]- sumW1;
        for(i = 1; i < kpoints; i++)
        {
            z[i] = z1[i-1];
            w[i] = w1[i-1];
        }
        
    }