Nektar++
Functions
Polylib Namespace Reference

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

Functions

double optdiff (double xl, double xr)
 The following function is used to circumvent/reduce "Subtractive Cancellation" The expression 1/dz is replaced by optinvsub(.,.) Added on 26 April 2017. More...
 
double laginterp (double z, int j, const double *zj, int np)
 
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. \( P_n^{\alpha,\beta}(z) = 0 \). 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 values and halves. 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 b of the orthogonal polynomials. More...
 
void JKMatrix (int n, double *a, double *b)
 Calcualtes the Jacobi-kronrod matrix by determining the a and coefficients. 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 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)\). 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 with the Gauss-Jacobi zeros. 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 with the Gauss-Radau-Jacobi zeros with a zero at z=-1. More...
 
void 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. More...
 
void 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. 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 the np Gauss-Jacobi points zgj at the arbitrary location z. 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 np Gauss-Radau-Jacobi points zgrj at the arbitrary location z. This routine assumes zgrj includes the point -1. 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 np Gauss-Radau-Jacobi points zgrj at the arbitrary location z. This routine assumes zgrj includes the point +1. 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 np Gauss-Lobatto-Jacobi points zgrj at the arbitrary location z. 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 arbitrary distribution at points zm. 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 (including z=-1) to an arbitrary distrubtion at points zm. 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 (including z=1) to an arbitrary distrubtion at points zm. 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 to an arbitrary distrubtion at points zm. 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 matrix from the three term recurrence relationship. More...
 
std::complex< Nektar::NekDouble > ImagBesselComp (int n, std::complex< Nektar::NekDouble > y)
 Calcualte the bessel function of the first kind with complex double input y. Taken from Numerical Recipies in C. More...
 

Detailed Description

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

Function Documentation

◆ chri1()

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 1522 of file Polylib.cpp.

Referenced by zwlk(), and zwrk().

1524  {
1525 
1526  double q = ceil(3.0*n/2);
1527  int size = (int)q+1;
1528  if(size < n+1)
1529  {
1530  fprintf(stderr,"input arrays a and b are too short\n");
1531  }
1532  double* r = new double[n+1];
1533  r[0] = z - a[0];
1534  r[1] = z - a[1] - b[1]/r[0];
1535  a0[0] = a[1] + r[1] - r[0];
1536  b0[0] = -r[0]*b[0];
1537 
1538  if(n == 1)
1539  {
1540  delete[] r;
1541  return;
1542  }
1543  int k = 0;
1544  for(k = 1; k < n; k++)
1545  {
1546  r[k+1] = z - a[k+1] - b[k+1]/r[k];
1547  a0[k] = a[k+1] + r[k+1] - r[k];
1548  b0[k] = b[k] * r[k]/r[k-1];
1549  }
1550  delete[] r;
1551 
1552  }

◆ Dgj()

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 543 of file Polylib.cpp.

References jacobd().

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

545  {
546 
547  double one = 1.0, two = 2.0;
548 
549  if (np <= 0){
550  D[0] = 0.0;
551  }
552  else{
553  int i,j;
554  double *pd;
555 
556  pd = (double *)malloc(np*sizeof(double));
557  jacobd(np,z,pd,np,alpha,beta);
558 
559  for (i = 0; i < np; i++){
560  for (j = 0; j < np; j++){
561 
562  if (i != j)
563  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
564  else
565  D[i*np+j] = (alpha - beta + (alpha + beta + two)*z[j])/
566  (two*(one - z[j]*z[j]));
567  }
568  }
569  free(pd);
570  }
571  return;
572  }
Polylib::jacobd
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.
Definition: Polylib.cpp:1131

◆ Dglj()

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 687 of file Polylib.cpp.

References gammaF(), and jacobd().

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

689  {
690  if (np <= 1){
691  D[0] = 0.0;
692  }
693  else{
694  int i, j;
695  double one = 1.0, two = 2.0;
696  double *pd;
697 
698  pd = (double *)malloc(np*sizeof(double));
699 
700  pd[0] = two*pow(-one,np)*gammaF(np + beta);
701  pd[0] /= gammaF(np - one)*gammaF(beta + two);
702  jacobd(np-2,z+1,pd+1,np-2,alpha+1,beta+1);
703  for(i = 1; i < np-1; ++i) pd[i] *= (one-z[i]*z[i]);
704  pd[np-1] = -two*gammaF(np + alpha);
705  pd[np-1] /= gammaF(np - one)*gammaF(alpha + two);
706 
707  for (i = 0; i < np; i++) {
708  for (j = 0; j < np; j++){
709  if (i != j)
710  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
711  else {
712  if (j == 0)
713  D[i*np+j] = (alpha - (np-1)*(np + alpha + beta))/(two*(beta+ two));
714  else if (j == np-1)
715  D[i*np+j] =-(beta - (np-1)*(np + alpha + beta))/(two*(alpha+ two));
716  else
717  D[i*np+j] = (alpha - beta + (alpha + beta)*z[j])/
718  (two*(one - z[j]*z[j]));
719  }
720  }
721  }
722  free(pd);
723  }
724 
725  return;
726  }
Polylib::jacobd
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.
Definition: Polylib.cpp:1131
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer values and halves.
Definition: Polylib.cpp:1158

◆ Dgrjm()

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 586 of file Polylib.cpp.

References gammaF(), and jacobd().

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

588  {
589 
590  if (np <= 0){
591  D[0] = 0.0;
592  }
593  else{
594  int i, j;
595  double one = 1.0, two = 2.0;
596  double *pd;
597 
598  pd = (double *)malloc(np*sizeof(double));
599 
600  pd[0] = pow(-one,np-1)*gammaF(np+beta+one);
601  pd[0] /= gammaF(np)*gammaF(beta+two);
602  jacobd(np-1,z+1,pd+1,np-1,alpha,beta+1);
603  for(i = 1; i < np; ++i) pd[i] *= (1+z[i]);
604 
605  for (i = 0; i < np; i++) {
606  for (j = 0; j < np; j++){
607  if (i != j)
608  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
609  else {
610  if(j == 0)
611  D[i*np+j] = -(np + alpha + beta + one)*(np - one)/
612  (two*(beta + two));
613  else
614  D[i*np+j] = (alpha - beta + one + (alpha + beta + one)*z[j])/
615  (two*(one - z[j]*z[j]));
616  }
617  }
618  }
619  free(pd);
620  }
621 
622  return;
623  }
Polylib::jacobd
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.
Definition: Polylib.cpp:1131
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer values and halves.
Definition: Polylib.cpp:1158

◆ Dgrjp()

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 636 of file Polylib.cpp.

References gammaF(), and jacobd().

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

638  {
639 
640  if (np <= 0){
641  D[0] = 0.0;
642  }
643  else{
644  int i, j;
645  double one = 1.0, two = 2.0;
646  double *pd;
647 
648  pd = (double *)malloc(np*sizeof(double));
649 
650 
651  jacobd(np-1,z,pd,np-1,alpha+1,beta);
652  for(i = 0; i < np-1; ++i) pd[i] *= (1-z[i]);
653  pd[np-1] = -gammaF(np+alpha+one);
654  pd[np-1] /= gammaF(np)*gammaF(alpha+two);
655 
656  for (i = 0; i < np; i++) {
657  for (j = 0; j < np; j++){
658  if (i != j)
659  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
660  else {
661  if(j == np-1)
662  D[i*np+j] = (np + alpha + beta + one)*(np - one)/
663  (two*(alpha + two));
664  else
665  D[i*np+j] = (alpha - beta - one + (alpha + beta + one)*z[j])/
666  (two*(one - z[j]*z[j]));
667  }
668  }
669  }
670  free(pd);
671  }
672 
673  return;
674  }
Polylib::jacobd
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.
Definition: Polylib.cpp:1131
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer values and halves.
Definition: Polylib.cpp:1158

◆ gammaF()

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 1158 of file Polylib.cpp.

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

1158  {
1159  double gamma = 1.0;
1160 
1161  if (x == -0.5) gamma = -2.0*sqrt(M_PI);
1162  else if (!x) return gamma;
1163  else if ((x-(int)x) == 0.5){
1164  int n = (int) x;
1165  double tmp = x;
1166 
1167  gamma = sqrt(M_PI);
1168  while(n--){
1169  tmp -= 1.0;
1170  gamma *= tmp;
1171  }
1172  }
1173  else if ((x-(int)x) == 0.0){
1174  int n = (int) x;
1175  double tmp = x;
1176 
1177  while(--n){
1178  tmp -= 1.0;
1179  gamma *= tmp;
1180  }
1181  }
1182  else
1183  fprintf(stderr,"%lf is not of integer or half order\n",x);
1184  return gamma;
1185  }

◆ hgj()

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 749 of file Polylib.cpp.

References EPS, and laginterp().

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

751  {
752  boost::ignore_unused(alpha, beta);
753  double zi, dz;
754 
755  zi = *(zgj+i);
756  dz = z-zi;
757  if (fabs(dz) < EPS) return 1.0;
758 
759  return laginterp(z, i, zgj, np);
760 
761  }
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:14
Polylib::laginterp
double laginterp(double z, int j, const double *zj, int np)
Definition: Polylib.cpp:43

◆ hglj()

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 855 of file Polylib.cpp.

References EPS, and laginterp().

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

857  {
858  boost::ignore_unused(alpha, beta);
859 
860  double zi, dz;
861 
862  zi = *(zglj+i);
863  dz = z-zi;
864  if (fabs(dz) < EPS) return 1.0;
865 
866  return laginterp(z, i, zglj, np);
867 
868  }
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:14
Polylib::laginterp
double laginterp(double z, int j, const double *zj, int np)
Definition: Polylib.cpp:43

◆ hgrjm()

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 784 of file Polylib.cpp.

References EPS, and laginterp().

Referenced by Imgrjm().

786  {
787  boost::ignore_unused(alpha, beta);
788 
789  double zi, dz;
790 
791  zi = *(zgrj+i);
792  dz = z-zi;
793  if (fabs(dz) < EPS) return 1.0;
794 
795  return laginterp(z, i, zgrj, np);
796  }
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:14
Polylib::laginterp
double laginterp(double z, int j, const double *zj, int np)
Definition: Polylib.cpp:43

◆ hgrjp()

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 820 of file Polylib.cpp.

References EPS, and laginterp().

Referenced by Imgrjp().

822  {
823  boost::ignore_unused(alpha, beta);
824  double zi, dz;
825 
826  zi = *(zgrj+i);
827  dz = z-zi;
828  if (fabs(dz) < EPS) return 1.0;
829 
830  return laginterp(z, i, zgrj, np);
831  }
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:14
Polylib::laginterp
double laginterp(double z, int j, const double *zj, int np)
Definition: Polylib.cpp:43

◆ ImagBesselComp()

std::complex< Nektar::NekDouble > Polylib::ImagBesselComp ( int  n,
std::complex< Nektar::NekDouble >  y 
)

Calcualte the bessel function of the first kind with complex double input y. Taken from Numerical Recipies in C.

Returns a complex double

Definition at line 1565 of file Polylib.cpp.

Referenced by Nektar::IncNavierStokes::SetUpWomersley().

1566  {
1567  std::complex<Nektar::NekDouble> z (1.0,0.0);
1568  std::complex<Nektar::NekDouble> zbes (1.0,0.0);
1569  std::complex<Nektar::NekDouble> zarg;
1570  Nektar::NekDouble tol = 1e-15;
1571  int maxit = 10000;
1572  int i = 1;
1573 
1574  zarg = -0.25*y*y;
1575 
1576  while (abs(z) > tol && i <= maxit){
1577  z = z*(1.0/i/(i+n)*zarg);
1578  if (abs(z) <= tol) break;
1579  zbes = zbes + z;
1580  i++;
1581  }
1582  zarg = 0.5*y;
1583  for (i=1;i<=n;i++){
1584  zbes = zbes*zarg;
1585  }
1586  return zbes;
1587 
1588  }
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43

◆ Imgj()

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 884 of file Polylib.cpp.

References hgj().

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

885  {
886  double zp;
887  int i, j;
888 
889  for (i = 0; i < nz; ++i) {
890  for (j = 0; j < mz; ++j)
891  {
892  zp = zm[j];
893  im [i*mz+j] = hgj(i, zp, zgj, nz, alpha, beta);
894  }
895  }
896 
897  return;
898  }
Polylib::hgj
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 the np Gauss-Jacobi points zgj at the ar...
Definition: Polylib.cpp:749

◆ Imglj()

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 974 of file Polylib.cpp.

References hglj().

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

976  {
977  double zp;
978  int i, j;
979 
980  for (i = 0; i < nz; i++) {
981  for (j = 0; j < mz; j++)
982  {
983  zp = zm[j];
984  im[i*mz+j] = hglj(i, zp, zglj, nz, alpha, beta);
985  }
986  }
987 
988  return;
989  }
Polylib::hglj
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 np Gauss-Lobatto-Jacobi points zgrj ...
Definition: Polylib.cpp:855

◆ Imgrjm()

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 913 of file Polylib.cpp.

References hgrjm().

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

915  {
916  double zp;
917  int i, j;
918 
919  for (i = 0; i < nz; i++) {
920  for (j = 0; j < mz; j++)
921  {
922  zp = zm[j];
923  im [i*mz+j] = hgrjm(i, zp, zgrj, nz, alpha, beta);
924  }
925  }
926 
927  return;
928  }
Polylib::hgrjm
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 np Gauss-Radau-Jacobi points zgrj at...
Definition: Polylib.cpp:784

◆ Imgrjp()

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 943 of file Polylib.cpp.

References hgrjp().

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

945  {
946  double zp;
947  int i, j;
948 
949  for (i = 0; i < nz; i++) {
950  for (j = 0; j < mz; j++)
951  {
952  zp = zm[j];
953  im [i*mz+j] = hgrjp(i, zp, zgrj, nz, alpha, beta);
954  }
955  }
956 
957  return;
958  }
Polylib::hgrjp
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 np Gauss-Radau-Jacobi points zgrj at...
Definition: Polylib.cpp:820

◆ jacobd()

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 1131 of file Polylib.cpp.

References jacobfd().

Referenced by Dgj(), Dglj(), Dgrjm(), Dgrjp(), Nektar::LibUtilities::NodalUtilTriangle::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilTetrahedron::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilPrism::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilQuad::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilHex::v_OrthoBasisDeriv(), Nektar::SolverUtils::AdvectionFR::v_SetupCFunctions(), Nektar::SolverUtils::DiffusionLFR::v_SetupCFunctions(), Nektar::SolverUtils::DiffusionLFRNS::v_SetupCFunctions(), and zwgj().

1133  {
1134  int i;
1135  double one = 1.0;
1136  if(n == 0)
1137  for(i = 0; i < np; ++i) polyd[i] = 0.0;
1138  else{
1139  //jacobf(np,z,polyd,n-1,alpha+one,beta+one);
1140  jacobfd(np,z,polyd,NULL,n-1,alpha+one,beta+one);
1141  for(i = 0; i < np; ++i) polyd[i] *= 0.5*(alpha + beta + (double)n + one);
1142  }
1143  return;
1144  }
Polylib::jacobfd
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, , and their first derivative, .
Definition: Polylib.cpp:1031

◆ jacobfd()

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 1031 of file Polylib.cpp.

Referenced by Nektar::LibUtilities::Basis::GenBasis(), jacobd(), Jacobz(), main(), Nektar::MultiRegions::ExpList1D::PeriodicEval(), Nektar::LibUtilities::NodalUtilTriangle::v_OrthoBasis(), Nektar::LibUtilities::NodalUtilTetrahedron::v_OrthoBasis(), Nektar::LibUtilities::NodalUtilPrism::v_OrthoBasis(), Nektar::LibUtilities::NodalUtilQuad::v_OrthoBasis(), Nektar::LibUtilities::NodalUtilHex::v_OrthoBasis(), Nektar::LibUtilities::NodalUtilTriangle::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilTetrahedron::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilPrism::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilQuad::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilHex::v_OrthoBasisDeriv(), zwglj(), zwgrjm(), and zwgrjp().

1032  {
1033  int i;
1034  double zero = 0.0, one = 1.0, two = 2.0;
1035 
1036  if(!np)
1037  return;
1038 
1039  if(n == 0){
1040  if(poly_in)
1041  for(i = 0; i < np; ++i)
1042  poly_in[i] = one;
1043  if(polyd)
1044  for(i = 0; i < np; ++i)
1045  polyd[i] = zero;
1046  }
1047  else if (n == 1){
1048  if(poly_in)
1049  for(i = 0; i < np; ++i)
1050  poly_in[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
1051  if(polyd)
1052  for(i = 0; i < np; ++i)
1053  polyd[i] = 0.5*(alpha + beta + two);
1054  }
1055  else{
1056  int k;
1057  double a1,a2,a3,a4;
1058  double two = 2.0, apb = alpha + beta;
1059  double *poly, *polyn1,*polyn2;
1060 
1061  if(poly_in){ // switch for case of no poynomial function return
1062  polyn1 = (double *)malloc(2*np*sizeof(double));
1063  polyn2 = polyn1+np;
1064  poly = poly_in;
1065  }
1066  else{
1067  polyn1 = (double *)malloc(3*np*sizeof(double));
1068  polyn2 = polyn1+np;
1069  poly = polyn2+np;
1070  }
1071 
1072  for(i = 0; i < np; ++i){
1073  polyn2[i] = one;
1074  polyn1[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
1075  }
1076 
1077  for(k = 2; k <= n; ++k){
1078  a1 = two*k*(k + apb)*(two*k + apb - two);
1079  a2 = (two*k + apb - one)*(alpha*alpha - beta*beta);
1080  a3 = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb);
1081  a4 = two*(k + alpha - one)*(k + beta - one)*(two*k + apb);
1082 
1083  a2 /= a1;
1084  a3 /= a1;
1085  a4 /= a1;
1086 
1087  for(i = 0; i < np; ++i){
1088  poly [i] = (a2 + a3*z[i])*polyn1[i] - a4*polyn2[i];
1089  polyn2[i] = polyn1[i];
1090  polyn1[i] = poly [i];
1091  }
1092  }
1093 
1094  if(polyd){
1095  a1 = n*(alpha - beta);
1096  a2 = n*(two*n + alpha + beta);
1097  a3 = two*(n + alpha)*(n + beta);
1098  a4 = (two*n + alpha + beta);
1099  a1 /= a4; a2 /= a4; a3 /= a4;
1100 
1101  // note polyn2 points to polyn1 at end of poly iterations
1102  for(i = 0; i < np; ++i){
1103  polyd[i] = (a1- a2*z[i])*poly[i] + a3*polyn2[i];
1104  polyd[i] /= (one - z[i]*z[i]);
1105  }
1106  }
1107 
1108  free(polyn1);
1109  }
1110 
1111  return;
1112  }

◆ Jacobz()

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 1195 of file Polylib.cpp.

References EPS, jacobfd(), and STOP.

1196  {
1197  int i,j,k;
1198  double dth = M_PI/(2.0*(double)n);
1199  double poly,pder,rlast=0.0;
1200  double sum,delr,r;
1201  double one = 1.0, two = 2.0;
1202 
1203  if(!n)
1204  return;
1205 
1206  for(k = 0; k < n; ++k){
1207  r = -cos((two*(double)k + one) * dth);
1208  if(k) r = 0.5*(r + rlast);
1209 
1210  for(j = 1; j < STOP; ++j){
1211  jacobfd(1,&r,&poly, &pder, n, alpha, beta);
1212 
1213  for(i = 0, sum = 0.0; i < k; ++i) sum += one/(r - z[i]);
1214 
1215  delr = -poly / (pder - sum * poly);
1216  r += delr;
1217  if( fabs(delr) < EPS ) break;
1218  }
1219  z[k] = r;
1220  rlast = r;
1221  }
1222  return;
1223  }
STOP
#define STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
Definition: Polylib.cpp:12
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:14
Polylib::jacobfd
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, , and their first derivative, .
Definition: Polylib.cpp:1031

◆ JacZeros()

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 1250 of file Polylib.cpp.

References RecCoeff(), and TriQL().

1251  {
1252 
1253  int i,j;
1254  RecCoeff(n,a,b,alpha,beta);
1255 
1256  double **z = new double*[n];
1257  for(i = 0; i < n; i++)
1258  {
1259  z[i] = new double[n];
1260  for(j = 0; j < n; j++)
1261  {
1262  z[i][j] = 0.0;
1263  }
1264  }
1265  for(i = 0; i < n; i++)
1266  {
1267  z[i][i] = 1.0;
1268  }
1269 
1270  // find eigenvalues and eigenvectors
1271  TriQL(n, a, b,z);
1272 
1273  delete[] z;
1274  return;
1275  }
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and b of the orthogonal polynomials.
Definition: Polylib.cpp:1281
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:1339

◆ JKMatrix()

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 1441 of file Polylib.cpp.

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

1442  {
1443  int i,j,k,m;
1444  // Working storage
1445  int size = (int)floor(n/2.0)+2;
1446  double *s = new double[size];
1447  double *t = new double[size];
1448 
1449  // Initialize s and t to zero
1450  for(i = 0; i < size; i++)
1451  {
1452  s[i] = 0.0;
1453  t[i] = 0.0;
1454  }
1455 
1456  t[1] = b[n+1];
1457  for(m = 0; m <= n-2; m++)
1458  {
1459  double u = 0.0;
1460  for(k = (int)floor((m+1)/2.0); k >= 0; k--)
1461  {
1462  int l = m-k;
1463  u = u+(a[k+n+1]-a[l])*t[k+1] + b[k+n+1]*s[k] - b[l]*s[k+1];
1464  s[k+1] = u;
1465  }
1466 
1467  // Swap the contents of s and t
1468  double *hold = s;
1469  s = t;
1470  t = hold;
1471  }
1472 
1473 
1474  for(j = (int)floor(n/2.0); j >= 0; j--)
1475  {
1476  s[j+1] = s[j];
1477  }
1478 
1479  for(m = n-1; m <= 2*n-3; m++)
1480  {
1481  double u = 0;
1482  for(k = m+1-n; k <= floor((m-1)/2.0); k++)
1483  {
1484  int l = m-k;
1485  j = n-1-l;
1486  u = u-(a[k+n+1]-a[l])*t[j+1] - b[k+n+1]*s[j+1] + b[l]*s[j+2];
1487  s[j+1] = u;
1488  }
1489 
1490  if(m%2 == 0)
1491  {
1492  k = m/2;
1493  a[k+n+1] = a[k] + (s[j+1]-b[k+n+1]*s[j+2])/t[j+2];
1494 
1495  }else
1496  {
1497  k = (m+1)/2;
1498  b[k+n+1] = s[j+1]/s[j+2];
1499  }
1500 
1501 
1502  // Swap the contents of s and t
1503  double *hold = s;
1504  s = t;
1505  t = hold;
1506  }
1507 
1508  a[2*n ] = a[n-1]-b[2*n]*s[1]/t[1];
1509 
1510  }

◆ laginterp()

double Polylib::laginterp ( double  z,
int  j,
const double *  zj,
int  np 
)

Definition at line 43 of file Polylib.cpp.

References optdiff().

Referenced by hgj(), hglj(), hgrjm(), and hgrjp().

44  {
45  double temp = 1.0;
46  for (int i=0; i<np; i++)
47  {
48  if (j != i)
49  {
50  temp *=optdiff(z,zj[i])/(zj[j]-zj[i]);
51  }
52  }
53  return temp;
54  }
Polylib::optdiff
double optdiff(double xl, double xr)
The following function is used to circumvent/reduce "Subtractive Cancellation" The expression 1/dz is...
Definition: Polylib.cpp:23

◆ optdiff()

double Polylib::optdiff ( double  xl,
double  xr 
)

The following function is used to circumvent/reduce "Subtractive Cancellation" The expression 1/dz is replaced by optinvsub(.,.) Added on 26 April 2017.

Definition at line 23 of file Polylib.cpp.

References L.

Referenced by laginterp().

24  {
25  double m_xln, m_xrn;
26  int m_expn;
27  int m_digits = static_cast<int>(fabs(floor(log10(DBL_EPSILON)))-1);
28 
29  if (fabs(xl-xr)<1.e-4){
30 
31  m_expn = static_cast<int>(floor(log10(fabs(xl-xr))));
32  m_xln = xl*powl(10.0L,-m_expn)-floor(xl*powl(10.0L,-m_expn)); // substract the digits overlap part
33  m_xrn = xr*powl(10.0L,-m_expn)-floor(xl*powl(10.0L,-m_expn)); // substract the common digits overlap part
34  m_xln = round(m_xln*powl(10.0L,m_digits+m_expn)); // git rid of rubbish
35  m_xrn = round(m_xrn*powl(10.0L,m_digits+m_expn));
36 
37  return powl(10.0L,-m_digits)*(m_xln-m_xrn);
38  }else{
39  return (xl-xr);
40  }
41  }
L
NekDouble L
Definition: CompressibleBL.cpp:79

◆ RecCoeff()

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 1281 of file Polylib.cpp.

References gammaF().

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

1282  {
1283 
1284  int i;
1285  double apb, apbi,a2b2;
1286 
1287  if(!n)
1288  return;
1289 
1290  // generate normalised terms
1291  apb = alpha + beta;
1292  apbi = 2.0 + apb;
1293 
1294  b[0] = pow(2.0,apb+1.0)*gammaF(alpha+1.0)*gammaF(beta+1.0)/gammaF(apbi); //MuZero
1295  a[0] = (beta-alpha)/apbi;
1296  b[1] = (4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));
1297 
1298  a2b2 = beta*beta-alpha*alpha;
1299 
1300  for(i = 1; i < n-1; i++){
1301  apbi = 2.0*(i+1) + apb;
1302  a[i] = a2b2/((apbi-2.0)*apbi);
1303  b[i+1] = (4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
1304  ((apbi*apbi-1)*apbi*apbi));
1305  }
1306 
1307  apbi = 2.0*n + apb;
1308  a[n-1] = a2b2/((apbi-2.0)*apbi);
1309 
1310  }
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer values and halves.
Definition: Polylib.cpp:1158

◆ TriQL()

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 1339 of file Polylib.cpp.

References CellMLToNektar.cellml_metadata::p, sign, and STOP.

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

1339  {
1340  int m,l,iter,i,k;
1341  double s,r,p,g,f,dd,c,b;
1342 
1343  double MuZero = e[0];
1344 
1345  // Renumber the elements of e
1346  for(i = 0; i < n-1; i++)
1347  {
1348  e[i] = sqrt(e[i+1]);
1349  }
1350  e[n-1] = 0.0;
1351 
1352 
1353  for (l=0;l<n;l++) {
1354  iter=0;
1355  do {
1356  for (m=l;m<n-1;m++) {
1357  dd=fabs(d[m])+fabs(d[m+1]);
1358  if (fabs(e[m])+dd == dd) break;
1359  }
1360  if (m != l) {
1361  if (iter++ == STOP){
1362  fprintf(stderr,"triQL: Too many iterations in TQLI");
1363  exit(1);
1364  }
1365  g=(d[l+1]-d[l])/(2.0*e[l]);
1366  r=sqrt((g*g)+1.0);
1367  g=d[m]-d[l]+e[l]/(g+sign(r,g));
1368  s=c=1.0;
1369  p=0.0;
1370  for (i=m-1;i>=l;i--) {
1371  f=s*e[i];
1372  b=c*e[i];
1373  if (fabs(f) >= fabs(g)) {
1374  c=g/f;
1375  r=sqrt((c*c)+1.0);
1376  e[i+1]=f*r;
1377  c *= (s=1.0/r);
1378  } else {
1379  s=f/g;
1380  r=sqrt((s*s)+1.0);
1381  e[i+1]=g*r;
1382  s *= (c=1.0/r);
1383  }
1384  g=d[i+1]-p;
1385  r=(d[i]-g)*s+2.0*c*b;
1386  p=s*r;
1387  d[i+1]=g+p;
1388  g=c*r-b;
1389 
1390  // Calculate the eigenvectors
1391  for(k = 0; k < n; k++)
1392  {
1393  f = z[k][i+1];
1394  z[k][i+1] = s*z[k][i] + c*f;
1395  z[k][i] = c*z[k][i] - s*f;
1396  }
1397 
1398  }
1399  d[l]=d[l]-p;
1400  e[l]=g;
1401  e[m]=0.0;
1402  }
1403  } while (m != l);
1404  }
1405 
1406  // order eigenvalues
1407  // Since we only need the first component of the eigenvectors
1408  // to calcualte the weight, we only swap the first components
1409  for(i = 0; i < n-1; ++i){
1410  k = i;
1411  p = d[i];
1412  for(l = i+1; l < n; ++l)
1413  if (d[l] < p) {
1414  k = l;
1415  p = d[l];
1416  }
1417  d[k] = d[i];
1418  d[i] = p;
1419 
1420  double temp = z[0][k];
1421  z[0][k] = z[0][i];
1422  z[0][i] = temp;
1423  }
1424 
1425  // Calculate the weights
1426  for(i =0 ; i < n; i++)
1427  {
1428  e[i] = MuZero*z[0][i]*z[0][i];
1429  }
1430  }
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:16
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
STOP
#define STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
Definition: Polylib.cpp:12

◆ zwgj()

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 91 of file Polylib.cpp.

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

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

93  {
94  int i;
95  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
96 
97  jacobz (np,z,alpha,beta);
98  jacobd (np,z,w,np,alpha,beta);
99 
100  fac = pow(two,apb + one)*gammaF(alpha + np + one)*gammaF(beta + np + one);
101  fac /= gammaF(np + one)*gammaF(apb + np + one);
102 
103  for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]*(one-z[i]*z[i]));
104 
105  return;
106  }
Polylib::jacobd
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.
Definition: Polylib.cpp:1131
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:60
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer values and halves.
Definition: Polylib.cpp:1158

◆ zwgk()

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 237 of file Polylib.cpp.

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

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

239  {
240 
241  int np = (npt-1)/2;
242 
243  int i,j;
244 
245  // number of kronrod points associated with the np gauss rule
246  int kpoints = 2*np + 1;
247 
248  // Define the number of required recurrence coefficents
249  int ncoeffs = (int)floor(3.0*(np+1)/2);
250 
251  // Define arrays for the recurrence coefficients
252  // We will use these arrays for the Kronrod results too, hence the
253  // reason for the size of the arrays
254  double *a = new double[kpoints];
255  double *b = new double[kpoints];
256 
257  // Initialize a and b to zero
258  for(i = 0; i < kpoints; i++)
259  {
260  a[i] = 0.0;
261  b[i] = 0.0;
262  }
263 
264  // Call the routine to calculate the recurrence coefficients
265  RecCoeff(ncoeffs,a,b,alpha,beta);
266 
267  // Call the routine to caluclate the jacobi-Kronrod matrix
268  JKMatrix(np,a,b);
269 
270  // Set up the identity matrix
271  double** zmatrix = new double*[kpoints];
272  for(i = 0; i < kpoints; i++)
273  {
274  zmatrix[i] = new double[kpoints];
275  for(j = 0; j < kpoints; j++)
276  {
277  zmatrix[i][j] = 0.0;
278  }
279  }
280  for(i = 0; i < kpoints; i++)
281  {
282  zmatrix[i][i] = 1.0;
283  }
284 
285  // Calculte the points and weights
286  TriQL(kpoints, a, b, zmatrix);
287 
288  for(i = 0; i < kpoints; i++)
289  {
290  z[i] = a[i];
291  w[i] = b[i];
292  }
293  delete[] a;
294  delete[] b;
295  for (i = 0; i < kpoints; i++)
296  {
297  delete[] zmatrix[i];
298  }
299  delete[] zmatrix;
300 
301  }
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and b of the orthogonal polynomials.
Definition: Polylib.cpp:1281
Polylib::JKMatrix
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the a and coefficients.
Definition: Polylib.cpp:1441
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:1339

◆ zwglj()

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 193 of file Polylib.cpp.

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

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

195  {
196 
197  if( np == 1 ){
198  z[0] = 0.0;
199  w[0] = 2.0;
200  }
201  else if( np == 2 ){
202  z[0] = -1.0;
203  z[1] = 1.0;
204 
205  w[0] = 1.0;
206  w[1] = 1.0;
207  }
208  else{
209  int i;
210  double fac, one = 1.0, apb = alpha + beta, two = 2.0;
211 
212  z[0] = -one;
213  z[np-1] = one;
214  jacobz (np-2,z + 1,alpha + one,beta + one);
215  jacobfd (np,z,w,NULL,np-1,alpha,beta);
216 
217  fac = pow(two,apb + 1)*gammaF(alpha + np)*gammaF(beta + np);
218  fac /= (np-1)*gammaF(np)*gammaF(alpha + beta + np + one);
219 
220  for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]);
221  w[0] *= (beta + one);
222  w[np-1] *= (alpha + one);
223  }
224 
225  return;
226  }
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:60
Polylib::jacobfd
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, , and their first derivative, .
Definition: Polylib.cpp:1031
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer values and halves.
Definition: Polylib.cpp:1158

◆ zwgrjm()

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 119 of file Polylib.cpp.

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

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

121  {
122 
123  if(np == 1){
124  z[0] = 0.0;
125  w[0] = 2.0;
126  }
127  else{
128  int i;
129  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
130 
131  z[0] = -one;
132  jacobz (np-1,z+1,alpha,beta+1);
133  jacobfd (np,z,w,NULL,np-1,alpha,beta);
134 
135  fac = pow(two,apb)*gammaF(alpha + np)*gammaF(beta + np);
136  fac /= gammaF(np)*(beta + np)*gammaF(apb + np + 1);
137 
138  for(i = 0; i < np; ++i) w[i] = fac*(1-z[i])/(w[i]*w[i]);
139  w[0] *= (beta + one);
140  }
141 
142  return;
143  }
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:60
Polylib::jacobfd
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, , and their first derivative, .
Definition: Polylib.cpp:1031
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer values and halves.
Definition: Polylib.cpp:1158

◆ zwgrjp()

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 157 of file Polylib.cpp.

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

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

159  {
160 
161  if(np == 1){
162  z[0] = 0.0;
163  w[0] = 2.0;
164  }
165  else{
166  int i;
167  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
168 
169  jacobz (np-1,z,alpha+1,beta);
170  z[np-1] = one;
171  jacobfd (np,z,w,NULL,np-1,alpha,beta);
172 
173  fac = pow(two,apb)*gammaF(alpha + np)*gammaF(beta + np);
174  fac /= gammaF(np)*(alpha + np)*gammaF(apb + np + 1);
175 
176  for(i = 0; i < np; ++i) w[i] = fac*(1+z[i])/(w[i]*w[i]);
177  w[np-1] *= (alpha + one);
178  }
179 
180  return;
181  }
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:60
Polylib::jacobfd
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, , and their first derivative, .
Definition: Polylib.cpp:1031
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer values and halves.
Definition: Polylib.cpp:1158

◆ zwlk()

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 417 of file Polylib.cpp.

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

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

419  {
420 
421  int np = (npt+1)/2;
422 
423  if(np < 4)
424  {
425  fprintf (stderr,"too few points in formula\n");
426  return;
427  }
428 
429  double endl = -1;
430  double endr = 1;
431  int i,j;
432 
433  // number of kronrod points associated with the np gauss rule
434  int kpoints = 2*np-1;
435 
436  // Define the number of required recurrence coefficents
437  int ncoeffs = (int)ceil(3.0*np/2)-1;
438 
439  // Define arrays for the recurrence coefficients
440  double *a = new double[ncoeffs+1];
441  double *b = new double[ncoeffs+1];
442 
443  // Initialize a and b to zero
444  for(i = 0; i < ncoeffs+1; i++)
445  {
446  a[i] = 0.0;
447  b[i] = 0.0;
448  }
449 
450  // Call the routine to calculate the recurrence coefficients
451  RecCoeff(ncoeffs,a,b,alpha,beta);
452 
453 
454  double* a0 = new double[ncoeffs];
455  double* b0 = new double[ncoeffs];
456 
457  chri1(ncoeffs,a,b,a0,b0,endl);
458 
459  double* a1 = new double[ncoeffs-1];
460  double* b1 = new double[ncoeffs-1];
461 
462  chri1(ncoeffs-1,a0,b0,a1,b1,endr);
463 
464 
465  double s = b1[0]/fabs(b1[0]);
466  b1[0] = s*b1[0];
467 
468  // Finding the 2*np-1 gauss-kronrod points
469  double* z1 = new double[2*np-3];
470  double* w1 = new double[2*np-3];
471  for(i = 0; i < ncoeffs; i++)
472  {
473  z1[i] = a1[i];
474  w1[i] = b1[i];
475  }
476  JKMatrix(np-2,z1,w1);
477  // Set up the identity matrix
478  double** zmatrix = new double*[2*np-3];
479  for(i = 0; i < 2*np-3; i++)
480  {
481  zmatrix[i] = new double[2*np-3];
482  for(j = 0; j < 2*np-3; j++)
483  {
484  zmatrix[i][j] = 0.0;
485  }
486  }
487  for(i = 0; i < 2*np-3; i++)
488  {
489  zmatrix[i][i] = 1.0;
490  }
491 
492  // Calculate the points and weights
493  TriQL(2*np-3, z1, w1, zmatrix);
494 
495  double sumW1 = 0.0;
496  double sumW1Z1 = 0.0;
497  for(i = 0; i < 2*np-3; i++)
498  {
499  w1[i] = s*w1[i]/(z1[i]-endl)/(z1[i]-endr);
500  sumW1 += w1[i];
501  sumW1Z1 += z1[i]*w1[i];
502  }
503 
504  double c0 = b[0]-sumW1;
505  double c1 = a[0]*b[0]-sumW1Z1;
506 
507  z[0] = endl;
508  z[2*np-2] = endr;
509  w[0] = (c0*endr-c1)/(endr-endl);
510  w[2*np-2] = (c1-c0*endl)/(endr-endl);
511 
512  for(i = 1; i < kpoints-1; i++)
513  {
514  z[i] = z1[i-1];
515  w[i] = w1[i-1];
516  }
517  delete[] a;
518  delete[] b;
519  delete[] a0;
520  delete[] b0;
521  delete[] a1;
522  delete[] b1;
523  delete[] z1;
524  delete[] w1;
525  for(i = 0; i < 2*np-3; i++)
526  {
527  delete[] zmatrix[i];
528  }
529  delete[] zmatrix;
530  }
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and b of the orthogonal polynomials.
Definition: Polylib.cpp:1281
Polylib::chri1
void chri1(int, double *, double *, double *, double *, double)
Given a weight function through the first n+1 coefficients a and b of its orthogonal polynomials thi...
Definition: Polylib.cpp:1522
Polylib::JKMatrix
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the a and coefficients.
Definition: Polylib.cpp:1441
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:1339

◆ zwrk()

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 310 of file Polylib.cpp.

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

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

312  {
313 
314  int np = npt/2;
315 
316  if(np < 2)
317  {
318  fprintf(stderr,"too few points in formula\n");
319  return;
320  }
321 
322  double end0 = -1;
323 
324  int i,j;
325 
326  // number of kronrod points associated with the np gauss rule
327  int kpoints = 2*np;
328 
329  // Define the number of required recurrence coefficents
330  int ncoeffs = (int)ceil(3.0*np/2);
331 
332  // Define arrays for the recurrence coefficients
333  double *a = new double[ncoeffs+1];
334  double *b = new double[ncoeffs+1];
335 
336  // Initialize a and b to zero
337  for(i = 0; i < ncoeffs+1; i++)
338  {
339  a[i] = 0.0;
340  b[i] = 0.0;
341  }
342 
343  // Call the routine to calculate the recurrence coefficients
344  RecCoeff(ncoeffs,a,b,alpha,beta);
345 
346  double* a0 = new double[ncoeffs];
347  double* b0 = new double[ncoeffs];
348 
349  chri1(ncoeffs,a,b,a0,b0,end0);
350 
351  double s = b0[0]/fabs(b0[0]);
352  b0[0] = s*b0[0];
353 
354  // Finding the 2*np-1 gauss-kronrod points
355  double* z1 = new double[2*np-1];
356  double* w1 = new double[2*np-1];
357  for(i = 0; i < ncoeffs; i++)
358  {
359  z1[i] = a0[i];
360  w1[i] = b0[i];
361  }
362  JKMatrix(np-1,z1,w1);
363  // Set up the identity matrix
364  double** zmatrix = new double*[2*np-1];
365  for(i = 0; i < 2*np-1; i++)
366  {
367  zmatrix[i] = new double[2*np-1];
368  for(j = 0; j < 2*np-1; j++)
369  {
370  zmatrix[i][j] = 0.0;
371  }
372  }
373  for(i = 0; i < 2*np-1; i++)
374  {
375  zmatrix[i][i] = 1.0;
376  }
377 
378  // Calculate the points and weights
379  TriQL(2*np-1, z1, w1, zmatrix);
380 
381  double sumW1 = 0.0;
382  for(i = 0; i < 2*np-1; i++)
383  {
384  w1[i] = s*w1[i]/(z1[i]-end0);
385  sumW1 += w1[i];
386  }
387 
388  z[0] = end0;
389  w[0] = b[0]- sumW1;
390  for(i = 1; i < kpoints; i++)
391  {
392  z[i] = z1[i-1];
393  w[i] = w1[i-1];
394  }
395 
396 
397  delete[] a;
398  delete[] b;
399  delete[] a0;
400  delete[] b0;
401  delete[] z1;
402  delete[] w1;
403  for(i = 0; i < 2*np-1; i++)
404  {
405  delete[] zmatrix[i];
406  }
407  delete[] zmatrix;
408  }
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and b of the orthogonal polynomials.
Definition: Polylib.cpp:1281
Polylib::chri1
void chri1(int, double *, double *, double *, double *, double)
Given a weight function through the first n+1 coefficients a and b of its orthogonal polynomials thi...
Definition: Polylib.cpp:1522
Polylib::JKMatrix
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the a and coefficients.
Definition: Polylib.cpp:1441
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:1339
Generated on Wed Jan 20 2021 03:34:22 for Nektar++ by   doxygen 1.8.13