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...
 
double gammaFracGammaF (const int x, const double alpha, const int y, const double beta)
 Calculate fraction of two Gamma functions, \( \Gamma(x+\alpha)/\Gamma(y+\beta) \), 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 1608 of file Polylib.cpp.

1610  {
1611 
1612  double q = ceil(3.0*n/2);
1613  int size = (int)q+1;
1614  if(size < n+1)
1615  {
1616  fprintf(stderr,"input arrays a and b are too short\n");
1617  }
1618  double* r = new double[n+1];
1619  r[0] = z - a[0];
1620  r[1] = z - a[1] - b[1]/r[0];
1621  a0[0] = a[1] + r[1] - r[0];
1622  b0[0] = -r[0]*b[0];
1623 
1624  if(n == 1)
1625  {
1626  delete[] r;
1627  return;
1628  }
1629  int k = 0;
1630  for(k = 1; k < n; k++)
1631  {
1632  r[k+1] = z - a[k+1] - b[k+1]/r[k];
1633  a0[k] = a[k+1] + r[k+1] - r[k];
1634  b0[k] = b[k] * r[k]/r[k-1];
1635  }
1636  delete[] r;
1637 
1638  }

Referenced by zwlk(), and zwrk().

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

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

References jacobd().

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

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

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

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

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

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

591  {
592 
593  if (np <= 0){
594  D[0] = 0.0;
595  }
596  else{
597  int i, j;
598  double one = 1.0, two = 2.0;
599  double *pd;
600 
601  pd = (double *)malloc(np*sizeof(double));
602 
603  pd[0] = pow(-one,np-1) * gammaFracGammaF(np + 1, beta, np, 0.0);
604  pd[0] /= gammaF(beta+two);
605  jacobd(np-1,z+1,pd+1,np-1,alpha,beta+1);
606  for(i = 1; i < np; ++i) pd[i] *= (1+z[i]);
607 
608  for (i = 0; i < np; i++) {
609  for (j = 0; j < np; j++){
610  if (i != j)
611  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
612  else {
613  if(j == 0)
614  D[i*np+j] = -(np + alpha + beta + one)*(np - one)/
615  (two*(beta + two));
616  else
617  D[i*np+j] = (alpha - beta + one + (alpha + beta + one)*z[j])/
618  (two*(one - z[j]*z[j]));
619  }
620  }
621  }
622  free(pd);
623  }
624 
625  return;
626  }

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

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

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

641  {
642 
643  if (np <= 0){
644  D[0] = 0.0;
645  }
646  else{
647  int i, j;
648  double one = 1.0, two = 2.0;
649  double *pd;
650 
651  pd = (double *)malloc(np*sizeof(double));
652 
653 
654  jacobd(np-1,z,pd,np-1,alpha+1,beta);
655  for(i = 0; i < np-1; ++i) pd[i] *= (1-z[i]);
656  pd[np-1] = -gammaFracGammaF(np + 1, alpha, np, 0.0);
657  pd[np-1] /= gammaF(alpha+two);
658 
659  for (i = 0; i < np; i++) {
660  for (j = 0; j < np; j++){
661  if (i != j)
662  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
663  else {
664  if(j == np-1)
665  D[i*np+j] = (np + alpha + beta + one)*(np - one)/
666  (two*(alpha + two));
667  else
668  D[i*np+j] = (alpha - beta - one + (alpha + beta + one)*z[j])/
669  (two*(one - z[j]*z[j]));
670  }
671  }
672  }
673  free(pd);
674  }
675 
676  return;
677  }

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

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

◆ gammaF()

double Polylib::gammaF ( const double  )

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

1161  {
1162  double gamma = 1.0;
1163 
1164  if (x == -0.5) gamma = -2.0*sqrt(M_PI);
1165  else if (!x) return gamma;
1166  else if ((x-(int)x) == 0.5){
1167  int n = (int) x;
1168  double tmp = x;
1169 
1170  gamma = sqrt(M_PI);
1171  while(n--){
1172  tmp -= 1.0;
1173  gamma *= tmp;
1174  }
1175  }
1176  else if ((x-(int)x) == 0.0){
1177  int n = (int) x;
1178  double tmp = x;
1179 
1180  while(--n){
1181  tmp -= 1.0;
1182  gamma *= tmp;
1183  }
1184  }
1185  else
1186  fprintf(stderr,"%lf is not of integer or half order\n",x);
1187  return gamma;
1188  }
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267

References tinysimd::sqrt().

Referenced by Dglj(), Dgrjm(), Dgrjp(), RecCoeff(), and test_gamma_fraction().

◆ gammaFracGammaF()

double Polylib::gammaFracGammaF ( const int  ,
const double  ,
const int  ,
const double   
)

Calculate fraction of two Gamma functions, \( \Gamma(x+\alpha)/\Gamma(y+\beta) \), 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)\)

Attempts simplification in cases like \( \Gamma(x+1)/\Gamma(x-1) = x*(x-1) \)

Definition at line 1203 of file Polylib.cpp.

1203  {
1204  double gamma = 1.0;
1205  double halfa = fabs(alpha - int(alpha));
1206  double halfb = fabs(beta - int(beta ));
1207  if (halfa == 0.0 && halfb == 0.0)// integer value
1208  {
1209  int X = x+alpha;
1210  int Y = y+beta;
1211  if (X > Y)
1212  {
1213  for(int tmp=X-1; tmp > Y-1; tmp-=1)
1214  gamma *= tmp;
1215  }
1216  else if (Y > X)
1217  {
1218  for(int tmp=Y-1; tmp > X-1; tmp-=1)
1219  gamma *= tmp;
1220  gamma = 1. / gamma;
1221  }
1222  }
1223  else if (halfa == 0.5 && halfb == 0.5)// both are halves
1224  {
1225  double X = x + alpha;
1226  double Y = y + beta;
1227  if( X > Y )
1228  {
1229  for(int tmp=int(X); tmp > int(Y); tmp-=1)
1230  gamma *= tmp-0.5;
1231  }
1232  else if( Y > X)
1233  {
1234  for(int tmp=int(Y); tmp > int(X); tmp-=1)
1235  gamma *= tmp-0.5;
1236  gamma = 1. / gamma;
1237  }
1238  }
1239  else
1240  {
1241  double X = x + alpha;
1242  double Y = y + beta;
1243  while (X>1 || Y>1)
1244  {
1245  if (X>1)
1246  {
1247  gamma *= X-1.;
1248  X -= 1.;
1249  }
1250  if (Y>1)
1251  {
1252  gamma /= Y-1.;
1253  Y -= 1.;
1254  }
1255  }
1256  if (X==0.5)
1257  {
1258  gamma *= sqrt(M_PI);
1259  }
1260  else if (Y==0.5)
1261  {
1262  gamma /= sqrt(M_PI);
1263  }
1264  else
1265  {
1266  fprintf(stderr,"%lf or %lf is not of integer or half order\n",X, Y);
1267  }
1268  }
1269 
1270  return gamma;
1271  }

References tinysimd::sqrt().

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

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

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

References EPS, and laginterp().

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

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

860  {
861  boost::ignore_unused(alpha, beta);
862 
863  double zi, dz;
864 
865  zi = *(zglj+i);
866  dz = z-zi;
867  if (fabs(dz) < EPS) return 1.0;
868 
869  return laginterp(z, i, zglj, np);
870 
871  }

References EPS, and laginterp().

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

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

789  {
790  boost::ignore_unused(alpha, beta);
791 
792  double zi, dz;
793 
794  zi = *(zgrj+i);
795  dz = z-zi;
796  if (fabs(dz) < EPS) return 1.0;
797 
798  return laginterp(z, i, zgrj, np);
799  }

References EPS, and laginterp().

Referenced by Imgrjm().

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

825  {
826  boost::ignore_unused(alpha, beta);
827  double zi, dz;
828 
829  zi = *(zgrj+i);
830  dz = z-zi;
831  if (fabs(dz) < EPS) return 1.0;
832 
833  return laginterp(z, i, zgrj, np);
834  }

References EPS, and laginterp().

Referenced by Imgrjp().

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

1652  {
1653  std::complex<Nektar::NekDouble> z (1.0,0.0);
1654  std::complex<Nektar::NekDouble> zbes (1.0,0.0);
1655  std::complex<Nektar::NekDouble> zarg;
1656  Nektar::NekDouble tol = 1e-15;
1657  int maxit = 10000;
1658  int i = 1;
1659 
1660  zarg = -0.25*y*y;
1661 
1662  while (abs(z) > tol && i <= maxit){
1663  z = z*(1.0/i/(i+n)*zarg);
1664  if (abs(z) <= tol) break;
1665  zbes = zbes + z;
1666  i++;
1667  }
1668  zarg = 0.5*y;
1669  for (i=1;i<=n;i++){
1670  zbes = zbes*zarg;
1671  }
1672  return zbes;
1673 
1674  }
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:272

References tinysimd::abs().

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

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

888  {
889  double zp;
890  int i, j;
891 
892  for (i = 0; i < nz; ++i) {
893  for (j = 0; j < mz; ++j)
894  {
895  zp = zm[j];
896  im [i*mz+j] = hgj(i, zp, zgj, nz, alpha, beta);
897  }
898  }
899 
900  return;
901  }
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:752

References hgj().

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

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

979  {
980  double zp;
981  int i, j;
982 
983  for (i = 0; i < nz; i++) {
984  for (j = 0; j < mz; j++)
985  {
986  zp = zm[j];
987  im[i*mz+j] = hglj(i, zp, zglj, nz, alpha, beta);
988  }
989  }
990 
991  return;
992  }
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:858

References hglj().

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

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

918  {
919  double zp;
920  int i, j;
921 
922  for (i = 0; i < nz; i++) {
923  for (j = 0; j < mz; j++)
924  {
925  zp = zm[j];
926  im [i*mz+j] = hgrjm(i, zp, zgrj, nz, alpha, beta);
927  }
928  }
929 
930  return;
931  }
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:787

References hgrjm().

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

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

948  {
949  double zp;
950  int i, j;
951 
952  for (i = 0; i < nz; i++) {
953  for (j = 0; j < mz; j++)
954  {
955  zp = zm[j];
956  im [i*mz+j] = hgrjp(i, zp, zgrj, nz, alpha, beta);
957  }
958  }
959 
960  return;
961  }
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:823

References hgrjp().

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

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

1136  {
1137  int i;
1138  double one = 1.0;
1139  if(n == 0)
1140  for(i = 0; i < np; ++i) polyd[i] = 0.0;
1141  else{
1142  //jacobf(np,z,polyd,n-1,alpha+one,beta+one);
1143  jacobfd(np,z,polyd,NULL,n-1,alpha+one,beta+one);
1144  for(i = 0; i < np; ++i) polyd[i] *= 0.5*(alpha + beta + (double)n + one);
1145  }
1146  return;
1147  }
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:1034

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().

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

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

Referenced by Nektar::LibUtilities::Basis::GenBasis(), jacobd(), Jacobz(), main(), 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().

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

1282  {
1283  int i,j,k;
1284  double dth = M_PI/(2.0*(double)n);
1285  double poly,pder,rlast=0.0;
1286  double sum,delr,r;
1287  double one = 1.0, two = 2.0;
1288 
1289  if(!n)
1290  return;
1291 
1292  for(k = 0; k < n; ++k){
1293  r = -cos((two*(double)k + one) * dth);
1294  if(k) r = 0.5*(r + rlast);
1295 
1296  for(j = 1; j < STOP; ++j){
1297  jacobfd(1,&r,&poly, &pder, n, alpha, beta);
1298 
1299  for(i = 0, sum = 0.0; i < k; ++i) sum += one/(r - z[i]);
1300 
1301  delr = -poly / (pder - sum * poly);
1302  r += delr;
1303  if( fabs(delr) < EPS ) break;
1304  }
1305  z[k] = r;
1306  rlast = r;
1307  }
1308  return;
1309  }
STOP
#define STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
Definition: Polylib.cpp:11

References EPS, jacobfd(), and STOP.

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

1337  {
1338 
1339  int i,j;
1340  RecCoeff(n,a,b,alpha,beta);
1341 
1342  double **z = new double*[n];
1343  for(i = 0; i < n; i++)
1344  {
1345  z[i] = new double[n];
1346  for(j = 0; j < n; j++)
1347  {
1348  z[i][j] = 0.0;
1349  }
1350  }
1351  for(i = 0; i < n; i++)
1352  {
1353  z[i][i] = 1.0;
1354  }
1355 
1356  // find eigenvalues and eigenvectors
1357  TriQL(n, a, b,z);
1358 
1359  delete[] z;
1360  return;
1361  }
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:1367
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:1425

References RecCoeff(), and TriQL().

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

1528  {
1529  int i,j,k,m;
1530  // Working storage
1531  int size = (int)floor(n/2.0)+2;
1532  double *s = new double[size];
1533  double *t = new double[size];
1534 
1535  // Initialize s and t to zero
1536  for(i = 0; i < size; i++)
1537  {
1538  s[i] = 0.0;
1539  t[i] = 0.0;
1540  }
1541 
1542  t[1] = b[n+1];
1543  for(m = 0; m <= n-2; m++)
1544  {
1545  double u = 0.0;
1546  for(k = (int)floor((m+1)/2.0); k >= 0; k--)
1547  {
1548  int l = m-k;
1549  u = u+(a[k+n+1]-a[l])*t[k+1] + b[k+n+1]*s[k] - b[l]*s[k+1];
1550  s[k+1] = u;
1551  }
1552 
1553  // Swap the contents of s and t
1554  double *hold = s;
1555  s = t;
1556  t = hold;
1557  }
1558 
1559 
1560  for(j = (int)floor(n/2.0); j >= 0; j--)
1561  {
1562  s[j+1] = s[j];
1563  }
1564 
1565  for(m = n-1; m <= 2*n-3; m++)
1566  {
1567  double u = 0;
1568  for(k = m+1-n; k <= floor((m-1)/2.0); k++)
1569  {
1570  int l = m-k;
1571  j = n-1-l;
1572  u = u-(a[k+n+1]-a[l])*t[j+1] - b[k+n+1]*s[j+1] + b[l]*s[j+2];
1573  s[j+1] = u;
1574  }
1575 
1576  if(m%2 == 0)
1577  {
1578  k = m/2;
1579  a[k+n+1] = a[k] + (s[j+1]-b[k+n+1]*s[j+2])/t[j+2];
1580 
1581  }else
1582  {
1583  k = (m+1)/2;
1584  b[k+n+1] = s[j+1]/s[j+2];
1585  }
1586 
1587 
1588  // Swap the contents of s and t
1589  double *hold = s;
1590  s = t;
1591  t = hold;
1592  }
1593 
1594  a[2*n ] = a[n-1]-b[2*n]*s[1]/t[1];
1595 
1596  }

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

◆ laginterp()

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

Definition at line 42 of file Polylib.cpp.

43  {
44  double temp = 1.0;
45  for (int i=0; i<np; i++)
46  {
47  if (j != i)
48  {
49  temp *=optdiff(z,zj[i])/(zj[j]-zj[i]);
50  }
51  }
52  return temp;
53  }
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:22

References optdiff().

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

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

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

Referenced by laginterp().

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

1368  {
1369 
1370  int i;
1371  double apb, apbi,a2b2;
1372 
1373  if(!n)
1374  return;
1375 
1376  // generate normalised terms
1377  apb = alpha + beta;
1378  apbi = 2.0 + apb;
1379 
1380  b[0] = pow(2.0,apb+1.0)*gammaF(alpha+1.0)*gammaF(beta+1.0)/gammaF(apbi); //MuZero
1381  a[0] = (beta-alpha)/apbi;
1382  b[1] = (4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));
1383 
1384  a2b2 = beta*beta-alpha*alpha;
1385 
1386  for(i = 1; i < n-1; i++){
1387  apbi = 2.0*(i+1) + apb;
1388  a[i] = a2b2/((apbi-2.0)*apbi);
1389  b[i+1] = (4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
1390  ((apbi*apbi-1)*apbi*apbi));
1391  }
1392 
1393  apbi = 2.0*n + apb;
1394  a[n-1] = a2b2/((apbi-2.0)*apbi);
1395 
1396  }

References gammaF().

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

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

1425  {
1426  int m,l,iter,i,k;
1427  double s,r,p,g,f,dd,c,b;
1428 
1429  double MuZero = e[0];
1430 
1431  // Renumber the elements of e
1432  for(i = 0; i < n-1; i++)
1433  {
1434  e[i] = sqrt(e[i+1]);
1435  }
1436  e[n-1] = 0.0;
1437 
1438 
1439  for (l=0;l<n;l++) {
1440  iter=0;
1441  do {
1442  for (m=l;m<n-1;m++) {
1443  dd=fabs(d[m])+fabs(d[m+1]);
1444  if (fabs(e[m])+dd == dd) break;
1445  }
1446  if (m != l) {
1447  if (iter++ == STOP){
1448  fprintf(stderr,"triQL: Too many iterations in TQLI");
1449  exit(1);
1450  }
1451  g=(d[l+1]-d[l])/(2.0*e[l]);
1452  r=sqrt((g*g)+1.0);
1453  g=d[m]-d[l]+e[l]/(g+sign(r,g));
1454  s=c=1.0;
1455  p=0.0;
1456  for (i=m-1;i>=l;i--) {
1457  f=s*e[i];
1458  b=c*e[i];
1459  if (fabs(f) >= fabs(g)) {
1460  c=g/f;
1461  r=sqrt((c*c)+1.0);
1462  e[i+1]=f*r;
1463  c *= (s=1.0/r);
1464  } else {
1465  s=f/g;
1466  r=sqrt((s*s)+1.0);
1467  e[i+1]=g*r;
1468  s *= (c=1.0/r);
1469  }
1470  g=d[i+1]-p;
1471  r=(d[i]-g)*s+2.0*c*b;
1472  p=s*r;
1473  d[i+1]=g+p;
1474  g=c*r-b;
1475 
1476  // Calculate the eigenvectors
1477  for(k = 0; k < n; k++)
1478  {
1479  f = z[k][i+1];
1480  z[k][i+1] = s*z[k][i] + c*f;
1481  z[k][i] = c*z[k][i] - s*f;
1482  }
1483 
1484  }
1485  d[l]=d[l]-p;
1486  e[l]=g;
1487  e[m]=0.0;
1488  }
1489  } while (m != l);
1490  }
1491 
1492  // order eigenvalues
1493  // Since we only need the first component of the eigenvectors
1494  // to calcualte the weight, we only swap the first components
1495  for(i = 0; i < n-1; ++i){
1496  k = i;
1497  p = d[i];
1498  for(l = i+1; l < n; ++l)
1499  if (d[l] < p) {
1500  k = l;
1501  p = d[l];
1502  }
1503  d[k] = d[i];
1504  d[i] = p;
1505 
1506  double temp = z[0][k];
1507  z[0][k] = z[0][i];
1508  z[0][i] = temp;
1509  }
1510 
1511  // Calculate the weights
1512  for(i =0 ; i < n; i++)
1513  {
1514  e[i] = MuZero*z[0][i]*z[0][i];
1515  }
1516  }
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350

References CellMLToNektar.cellml_metadata::p, sign, tinysimd::sqrt(), and STOP.

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

◆ 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.

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) * gammaFracGammaF(np + 1, alpha, np + 1, 0.0)
101  * gammaFracGammaF(np + 1, beta , np + 1, apb);
102 
103  for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]*(one-z[i]*z[i]));
104 
105  return;
106  }
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:59

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

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

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

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

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

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

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

197  {
198 
199  if( np == 1 ){
200  z[0] = 0.0;
201  w[0] = 2.0;
202  }
203  else if( np == 2 ){
204  z[0] = -1.0;
205  z[1] = 1.0;
206 
207  w[0] = 1.0;
208  w[1] = 1.0;
209  }
210  else{
211  int i;
212  double fac, one = 1.0, apb = alpha + beta, two = 2.0;
213 
214  z[0] = -one;
215  z[np-1] = one;
216  jacobz (np-2,z + 1,alpha + one,beta + one);
217  jacobfd (np,z,w,NULL,np-1,alpha,beta);
218 
219  fac = pow(two,apb + 1) * gammaFracGammaF(np, alpha, np, 0.0)
220  * gammaFracGammaF(np, beta , np + 1, apb);
221  fac /= (np-1);
222 
223  for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]);
224  w[0] *= (beta + one);
225  w[np-1] *= (alpha + one);
226  }
227 
228  return;
229  }

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

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

◆ 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.

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) * gammaFracGammaF(np, alpha, np, 0.0)
136  * gammaFracGammaF(np, beta , np + 1, apb);
137  fac /= (beta + np);
138 
139  for(i = 0; i < np; ++i) w[i] = fac*(1-z[i])/(w[i]*w[i]);
140  w[0] *= (beta + one);
141  }
142 
143  return;
144  }

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

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

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

160  {
161 
162  if(np == 1){
163  z[0] = 0.0;
164  w[0] = 2.0;
165  }
166  else{
167  int i;
168  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
169 
170  jacobz (np-1,z,alpha+1,beta);
171  z[np-1] = one;
172  jacobfd (np,z,w,NULL,np-1,alpha,beta);
173 
174  fac = pow(two,apb) * gammaFracGammaF(np, alpha, np, 0.0)
175  * gammaFracGammaF(np, beta , np + 1, apb);
176  fac /= (alpha + np);
177 
178  for(i = 0; i < np; ++i) w[i] = fac*(1+z[i])/(w[i]*w[i]);
179  w[np-1] *= (alpha + one);
180  }
181 
182  return;
183  }

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

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

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

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

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

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

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

315  {
316 
317  int np = npt/2;
318 
319  if(np < 2)
320  {
321  fprintf(stderr,"too few points in formula\n");
322  return;
323  }
324 
325  double end0 = -1;
326 
327  int i,j;
328 
329  // number of kronrod points associated with the np gauss rule
330  int kpoints = 2*np;
331 
332  // Define the number of required recurrence coefficents
333  int ncoeffs = (int)ceil(3.0*np/2);
334 
335  // Define arrays for the recurrence coefficients
336  double *a = new double[ncoeffs+1];
337  double *b = new double[ncoeffs+1];
338 
339  // Initialize a and b to zero
340  for(i = 0; i < ncoeffs+1; i++)
341  {
342  a[i] = 0.0;
343  b[i] = 0.0;
344  }
345 
346  // Call the routine to calculate the recurrence coefficients
347  RecCoeff(ncoeffs,a,b,alpha,beta);
348 
349  double* a0 = new double[ncoeffs];
350  double* b0 = new double[ncoeffs];
351 
352  chri1(ncoeffs,a,b,a0,b0,end0);
353 
354  double s = b0[0]/fabs(b0[0]);
355  b0[0] = s*b0[0];
356 
357  // Finding the 2*np-1 gauss-kronrod points
358  double* z1 = new double[2*np-1];
359  double* w1 = new double[2*np-1];
360  for(i = 0; i < ncoeffs; i++)
361  {
362  z1[i] = a0[i];
363  w1[i] = b0[i];
364  }
365  JKMatrix(np-1,z1,w1);
366  // Set up the identity matrix
367  double** zmatrix = new double*[2*np-1];
368  for(i = 0; i < 2*np-1; i++)
369  {
370  zmatrix[i] = new double[2*np-1];
371  for(j = 0; j < 2*np-1; j++)
372  {
373  zmatrix[i][j] = 0.0;
374  }
375  }
376  for(i = 0; i < 2*np-1; i++)
377  {
378  zmatrix[i][i] = 1.0;
379  }
380 
381  // Calculate the points and weights
382  TriQL(2*np-1, z1, w1, zmatrix);
383 
384  double sumW1 = 0.0;
385  for(i = 0; i < 2*np-1; i++)
386  {
387  w1[i] = s*w1[i]/(z1[i]-end0);
388  sumW1 += w1[i];
389  }
390 
391  z[0] = end0;
392  w[0] = b[0]- sumW1;
393  for(i = 1; i < kpoints; i++)
394  {
395  z[i] = z1[i-1];
396  w[i] = w1[i-1];
397  }
398 
399 
400  delete[] a;
401  delete[] b;
402  delete[] a0;
403  delete[] b0;
404  delete[] z1;
405  delete[] w1;
406  for(i = 0; i < 2*np-1; i++)
407  {
408  delete[] zmatrix[i];
409  }
410  delete[] zmatrix;
411  }

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

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

Generated on Wed Nov 24 2021 10:33:38 for Nektar++ by doxygen 1.9.1