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

1692 {
1693 
1694  double q = ceil(3.0 * n / 2);
1695  int size = (int)q + 1;
1696  if (size < n + 1)
1697  {
1698  fprintf(stderr, "input arrays a and b are too short\n");
1699  }
1700  double *r = new double[n + 1];
1701  r[0] = z - a[0];
1702  r[1] = z - a[1] - b[1] / r[0];
1703  a0[0] = a[1] + r[1] - r[0];
1704  b0[0] = -r[0] * b[0];
1705 
1706  if (n == 1)
1707  {
1708  delete[] r;
1709  return;
1710  }
1711  int k = 0;
1712  for (k = 1; k < n; k++)
1713  {
1714  r[k + 1] = z - a[k + 1] - b[k + 1] / r[k];
1715  a0[k] = a[k + 1] + r[k + 1] - r[k];
1716  b0[k] = b[k] * r[k] / r[k - 1];
1717  }
1718  delete[] r;
1719 }

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

560 {
561 
562  double one = 1.0, two = 2.0;
563 
564  if (np <= 0)
565  {
566  D[0] = 0.0;
567  }
568  else
569  {
570  int i, j;
571  double *pd;
572 
573  pd = (double *)malloc(np * sizeof(double));
574  jacobd(np, z, pd, np, alpha, beta);
575 
576  for (i = 0; i < np; i++)
577  {
578  for (j = 0; j < np; j++)
579  {
580 
581  if (i != j)
582  D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
583  else
584  D[i * np + j] =
585  (alpha - beta + (alpha + beta + two) * z[j]) /
586  (two * (one - z[j] * z[j]));
587  }
588  }
589  free(pd);
590  }
591  return;
592 }
Nektar::LibUtilities::beta
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:61
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:1187

References Nektar::LibUtilities::beta, and 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 718 of file Polylib.cpp.

720 {
721  if (np <= 1)
722  {
723  D[0] = 0.0;
724  }
725  else
726  {
727  int i, j;
728  double one = 1.0, two = 2.0;
729  double *pd;
730 
731  pd = (double *)malloc(np * sizeof(double));
732 
733  pd[0] = two * pow(-one, np) * gammaFracGammaF(np, beta, np - 1, 0.0);
734  pd[0] /= gammaF(beta + two);
735  jacobd(np - 2, z + 1, pd + 1, np - 2, alpha + 1, beta + 1);
736  for (i = 1; i < np - 1; ++i)
737  pd[i] *= (one - z[i] * z[i]);
738  pd[np - 1] = -two * gammaFracGammaF(np, alpha, np - 1, 0.0);
739  pd[np - 1] /= gammaF(alpha + two);
740 
741  for (i = 0; i < np; i++)
742  {
743  for (j = 0; j < np; j++)
744  {
745  if (i != j)
746  D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
747  else
748  {
749  if (j == 0)
750  D[i * np + j] =
751  (alpha - (np - 1) * (np + alpha + beta)) /
752  (two * (beta + two));
753  else if (j == np - 1)
754  D[i * np + j] =
755  -(beta - (np - 1) * (np + alpha + beta)) /
756  (two * (alpha + two));
757  else
758  D[i * np + j] = (alpha - beta + (alpha + beta) * z[j]) /
759  (two * (one - z[j] * z[j]));
760  }
761  }
762  }
763  free(pd);
764  }
765 
766  return;
767 }
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer values and halves.
Definition: Polylib.cpp:1216
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:1265

References Nektar::LibUtilities::beta, 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 605 of file Polylib.cpp.

607 {
608 
609  if (np <= 0)
610  {
611  D[0] = 0.0;
612  }
613  else
614  {
615  int i, j;
616  double one = 1.0, two = 2.0;
617  double *pd;
618 
619  pd = (double *)malloc(np * sizeof(double));
620 
621  pd[0] = pow(-one, np - 1) * gammaFracGammaF(np + 1, beta, np, 0.0);
622  pd[0] /= gammaF(beta + two);
623  jacobd(np - 1, z + 1, pd + 1, np - 1, alpha, beta + 1);
624  for (i = 1; i < np; ++i)
625  pd[i] *= (1 + z[i]);
626 
627  for (i = 0; i < np; i++)
628  {
629  for (j = 0; j < np; j++)
630  {
631  if (i != j)
632  D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
633  else
634  {
635  if (j == 0)
636  D[i * np + j] = -(np + alpha + beta + one) *
637  (np - one) / (two * (beta + two));
638  else
639  D[i * np + j] =
640  (alpha - beta + one + (alpha + beta + one) * z[j]) /
641  (two * (one - z[j] * z[j]));
642  }
643  }
644  }
645  free(pd);
646  }
647 
648  return;
649 }

References Nektar::LibUtilities::beta, 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 661 of file Polylib.cpp.

663 {
664 
665  if (np <= 0)
666  {
667  D[0] = 0.0;
668  }
669  else
670  {
671  int i, j;
672  double one = 1.0, two = 2.0;
673  double *pd;
674 
675  pd = (double *)malloc(np * sizeof(double));
676 
677  jacobd(np - 1, z, pd, np - 1, alpha + 1, beta);
678  for (i = 0; i < np - 1; ++i)
679  pd[i] *= (1 - z[i]);
680  pd[np - 1] = -gammaFracGammaF(np + 1, alpha, np, 0.0);
681  pd[np - 1] /= gammaF(alpha + two);
682 
683  for (i = 0; i < np; i++)
684  {
685  for (j = 0; j < np; j++)
686  {
687  if (i != j)
688  D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
689  else
690  {
691  if (j == np - 1)
692  D[i * np + j] = (np + alpha + beta + one) * (np - one) /
693  (two * (alpha + two));
694  else
695  D[i * np + j] =
696  (alpha - beta - one + (alpha + beta + one) * z[j]) /
697  (two * (one - z[j] * z[j]));
698  }
699  }
700  }
701  free(pd);
702  }
703 
704  return;
705 }

References Nektar::LibUtilities::beta, 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 1216 of file Polylib.cpp.

1217 {
1218  double gamma = 1.0;
1219 
1220  if (x == -0.5)
1221  gamma = -2.0 * sqrt(M_PI);
1222  else if (!x)
1223  return gamma;
1224  else if ((x - (int)x) == 0.5)
1225  {
1226  int n = (int)x;
1227  double tmp = x;
1228 
1229  gamma = sqrt(M_PI);
1230  while (n--)
1231  {
1232  tmp -= 1.0;
1233  gamma *= tmp;
1234  }
1235  }
1236  else if ((x - (int)x) == 0.0)
1237  {
1238  int n = (int)x;
1239  double tmp = x;
1240 
1241  while (--n)
1242  {
1243  tmp -= 1.0;
1244  gamma *= tmp;
1245  }
1246  }
1247  else
1248  fprintf(stderr, "%lf is not of integer or half order\n", x);
1249  return gamma;
1250 }
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291

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

1267 {
1268  double gamma = 1.0;
1269  double halfa = fabs(alpha - int(alpha));
1270  double halfb = fabs(beta - int(beta));
1271  if (halfa == 0.0 && halfb == 0.0) // integer value
1272  {
1273  int X = x + alpha;
1274  int Y = y + beta;
1275  if (X > Y)
1276  {
1277  for (int tmp = X - 1; tmp > Y - 1; tmp -= 1)
1278  gamma *= tmp;
1279  }
1280  else if (Y > X)
1281  {
1282  for (int tmp = Y - 1; tmp > X - 1; tmp -= 1)
1283  gamma *= tmp;
1284  gamma = 1. / gamma;
1285  }
1286  }
1287  else if (halfa == 0.5 && halfb == 0.5) // both are halves
1288  {
1289  double X = x + alpha;
1290  double Y = y + beta;
1291  if (X > Y)
1292  {
1293  for (int tmp = int(X); tmp > int(Y); tmp -= 1)
1294  gamma *= tmp - 0.5;
1295  }
1296  else if (Y > X)
1297  {
1298  for (int tmp = int(Y); tmp > int(X); tmp -= 1)
1299  gamma *= tmp - 0.5;
1300  gamma = 1. / gamma;
1301  }
1302  }
1303  else
1304  {
1305  double X = x + alpha;
1306  double Y = y + beta;
1307  while (X > 1 || Y > 1)
1308  {
1309  if (X > 1)
1310  {
1311  gamma *= X - 1.;
1312  X -= 1.;
1313  }
1314  if (Y > 1)
1315  {
1316  gamma /= Y - 1.;
1317  Y -= 1.;
1318  }
1319  }
1320  if (X == 0.5)
1321  {
1322  gamma *= sqrt(M_PI);
1323  }
1324  else if (Y == 0.5)
1325  {
1326  gamma /= sqrt(M_PI);
1327  }
1328  else
1329  {
1330  fprintf(stderr, "%lf or %lf is not of integer or half order\n", X,
1331  Y);
1332  }
1333  }
1334 
1335  return gamma;
1336 }

References Nektar::LibUtilities::beta, and 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 789 of file Polylib.cpp.

791 {
792  boost::ignore_unused(alpha, beta);
793  double zi, dz;
794 
795  zi = *(zgj + i);
796  dz = z - zi;
797  if (fabs(dz) < EPS)
798  return 1.0;
799 
800  return laginterp(z, i, zgj, np);
801 }
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:53

References Nektar::LibUtilities::beta, 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 895 of file Polylib.cpp.

897 {
898  boost::ignore_unused(alpha, beta);
899 
900  double zi, dz;
901 
902  zi = *(zglj + i);
903  dz = z - zi;
904  if (fabs(dz) < EPS)
905  return 1.0;
906 
907  return laginterp(z, i, zglj, np);
908 }

References Nektar::LibUtilities::beta, 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 824 of file Polylib.cpp.

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

References Nektar::LibUtilities::beta, 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 860 of file Polylib.cpp.

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

References Nektar::LibUtilities::beta, 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 1731 of file Polylib.cpp.

1733 {
1734  std::complex<Nektar::NekDouble> z(1.0, 0.0);
1735  std::complex<Nektar::NekDouble> zbes(1.0, 0.0);
1736  std::complex<Nektar::NekDouble> zarg;
1737  Nektar::NekDouble tol = 1e-15;
1738  int maxit = 10000;
1739  int i = 1;
1740 
1741  zarg = -0.25 * y * y;
1742 
1743  while (abs(z) > tol && i <= maxit)
1744  {
1745  z = z * (1.0 / i / (i + n) * zarg);
1746  if (abs(z) <= tol)
1747  break;
1748  zbes = zbes + z;
1749  i++;
1750  }
1751  zarg = 0.5 * y;
1752  for (i = 1; i <= n; i++)
1753  {
1754  zbes = zbes * zarg;
1755  }
1756  return zbes;
1757 }
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:295

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

925 {
926  double zp;
927  int i, j;
928 
929  for (i = 0; i < nz; ++i)
930  {
931  for (j = 0; j < mz; ++j)
932  {
933  zp = zm[j];
934  im[i * mz + j] = hgj(i, zp, zgj, nz, alpha, beta);
935  }
936  }
937 
938  return;
939 }
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:789

References Nektar::LibUtilities::beta, and 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 1016 of file Polylib.cpp.

1018 {
1019  double zp;
1020  int i, j;
1021 
1022  for (i = 0; i < nz; i++)
1023  {
1024  for (j = 0; j < mz; j++)
1025  {
1026  zp = zm[j];
1027  im[i * mz + j] = hglj(i, zp, zglj, nz, alpha, beta);
1028  }
1029  }
1030 
1031  return;
1032 }
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:895

References Nektar::LibUtilities::beta, and 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 954 of file Polylib.cpp.

956 {
957  double zp;
958  int i, j;
959 
960  for (i = 0; i < nz; i++)
961  {
962  for (j = 0; j < mz; j++)
963  {
964  zp = zm[j];
965  im[i * mz + j] = hgrjm(i, zp, zgrj, nz, alpha, beta);
966  }
967  }
968 
969  return;
970 }
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:824

References Nektar::LibUtilities::beta, and 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 985 of file Polylib.cpp.

987 {
988  double zp;
989  int i, j;
990 
991  for (i = 0; i < nz; i++)
992  {
993  for (j = 0; j < mz; j++)
994  {
995  zp = zm[j];
996  im[i * mz + j] = hgrjp(i, zp, zgrj, nz, alpha, beta);
997  }
998  }
999 
1000  return;
1001 }
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:860

References Nektar::LibUtilities::beta, and 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 1187 of file Polylib.cpp.

1189 {
1190  int i;
1191  double one = 1.0;
1192  if (n == 0)
1193  for (i = 0; i < np; ++i)
1194  polyd[i] = 0.0;
1195  else
1196  {
1197  // jacobf(np,z,polyd,n-1,alpha+one,beta+one);
1198  jacobfd(np, z, polyd, NULL, n - 1, alpha + one, beta + one);
1199  for (i = 0; i < np; ++i)
1200  polyd[i] *= 0.5 * (alpha + beta + (double)n + one);
1201  }
1202  return;
1203 }
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:1074

References Nektar::LibUtilities::beta, and 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 1074 of file Polylib.cpp.

1076 {
1077  int i;
1078  double zero = 0.0, one = 1.0, two = 2.0;
1079 
1080  if (!np)
1081  return;
1082 
1083  if (n == 0)
1084  {
1085  if (poly_in)
1086  for (i = 0; i < np; ++i)
1087  poly_in[i] = one;
1088  if (polyd)
1089  for (i = 0; i < np; ++i)
1090  polyd[i] = zero;
1091  }
1092  else if (n == 1)
1093  {
1094  if (poly_in)
1095  for (i = 0; i < np; ++i)
1096  poly_in[i] = 0.5 * (alpha - beta + (alpha + beta + two) * z[i]);
1097  if (polyd)
1098  for (i = 0; i < np; ++i)
1099  polyd[i] = 0.5 * (alpha + beta + two);
1100  }
1101  else
1102  {
1103  int k;
1104  double a1, a2, a3, a4;
1105  double two = 2.0, apb = alpha + beta;
1106  double *poly, *polyn1, *polyn2;
1107 
1108  if (poly_in)
1109  { // switch for case of no poynomial function return
1110  polyn1 = (double *)malloc(2 * np * sizeof(double));
1111  polyn2 = polyn1 + np;
1112  poly = poly_in;
1113  }
1114  else
1115  {
1116  polyn1 = (double *)malloc(3 * np * sizeof(double));
1117  polyn2 = polyn1 + np;
1118  poly = polyn2 + np;
1119  }
1120 
1121  for (i = 0; i < np; ++i)
1122  {
1123  polyn2[i] = one;
1124  polyn1[i] = 0.5 * (alpha - beta + (alpha + beta + two) * z[i]);
1125  }
1126 
1127  for (k = 2; k <= n; ++k)
1128  {
1129  a1 = two * k * (k + apb) * (two * k + apb - two);
1130  a2 = (two * k + apb - one) * (alpha * alpha - beta * beta);
1131  a3 =
1132  (two * k + apb - two) * (two * k + apb - one) * (two * k + apb);
1133  a4 = two * (k + alpha - one) * (k + beta - one) * (two * k + apb);
1134 
1135  a2 /= a1;
1136  a3 /= a1;
1137  a4 /= a1;
1138 
1139  for (i = 0; i < np; ++i)
1140  {
1141  poly[i] = (a2 + a3 * z[i]) * polyn1[i] - a4 * polyn2[i];
1142  polyn2[i] = polyn1[i];
1143  polyn1[i] = poly[i];
1144  }
1145  }
1146 
1147  if (polyd)
1148  {
1149  a1 = n * (alpha - beta);
1150  a2 = n * (two * n + alpha + beta);
1151  a3 = two * (n + alpha) * (n + beta);
1152  a4 = (two * n + alpha + beta);
1153  a1 /= a4;
1154  a2 /= a4;
1155  a3 /= a4;
1156 
1157  // note polyn2 points to polyn1 at end of poly iterations
1158  for (i = 0; i < np; ++i)
1159  {
1160  polyd[i] = (a1 - a2 * z[i]) * poly[i] + a3 * polyn2[i];
1161  polyd[i] /= (one - z[i] * z[i]);
1162  }
1163  }
1164 
1165  free(polyn1);
1166  }
1167 
1168  return;
1169 }

References Nektar::LibUtilities::beta.

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

1348 {
1349  int i, j, k;
1350  double dth = M_PI / (2.0 * (double)n);
1351  double poly, pder, rlast = 0.0;
1352  double sum, delr, r;
1353  double one = 1.0, two = 2.0;
1354 
1355  if (!n)
1356  return;
1357 
1358  for (k = 0; k < n; ++k)
1359  {
1360  r = -cos((two * (double)k + one) * dth);
1361  if (k)
1362  r = 0.5 * (r + rlast);
1363 
1364  for (j = 1; j < STOP; ++j)
1365  {
1366  jacobfd(1, &r, &poly, &pder, n, alpha, beta);
1367 
1368  for (i = 0, sum = 0.0; i < k; ++i)
1369  sum += one / (r - z[i]);
1370 
1371  delr = -poly / (pder - sum * poly);
1372  r += delr;
1373  if (fabs(delr) < EPS)
1374  break;
1375  }
1376  z[k] = r;
1377  rlast = r;
1378  }
1379  return;
1380 }
STOP
#define STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
Definition: Polylib.cpp:11

References Nektar::LibUtilities::beta, 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 1406 of file Polylib.cpp.

1408 {
1409 
1410  int i, j;
1411  RecCoeff(n, a, b, alpha, beta);
1412 
1413  double **z = new double *[n];
1414  for (i = 0; i < n; i++)
1415  {
1416  z[i] = new double[n];
1417  for (j = 0; j < n; j++)
1418  {
1419  z[i][j] = 0.0;
1420  }
1421  }
1422  for (i = 0; i < n; i++)
1423  {
1424  z[i][i] = 1.0;
1425  }
1426 
1427  // find eigenvalues and eigenvectors
1428  TriQL(n, a, b, z);
1429 
1430  delete[] z;
1431  return;
1432 }
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:1438
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:1497

References Nektar::LibUtilities::beta, 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 1610 of file Polylib.cpp.

1611 {
1612  int i, j, k, m;
1613  // Working storage
1614  int size = (int)floor(n / 2.0) + 2;
1615  double *s = new double[size];
1616  double *t = new double[size];
1617 
1618  // Initialize s and t to zero
1619  for (i = 0; i < size; i++)
1620  {
1621  s[i] = 0.0;
1622  t[i] = 0.0;
1623  }
1624 
1625  t[1] = b[n + 1];
1626  for (m = 0; m <= n - 2; m++)
1627  {
1628  double u = 0.0;
1629  for (k = (int)floor((m + 1) / 2.0); k >= 0; k--)
1630  {
1631  int l = m - k;
1632  u = u + (a[k + n + 1] - a[l]) * t[k + 1] + b[k + n + 1] * s[k] -
1633  b[l] * s[k + 1];
1634  s[k + 1] = u;
1635  }
1636 
1637  // Swap the contents of s and t
1638  double *hold = s;
1639  s = t;
1640  t = hold;
1641  }
1642 
1643  for (j = (int)floor(n / 2.0); j >= 0; j--)
1644  {
1645  s[j + 1] = s[j];
1646  }
1647 
1648  for (m = n - 1; m <= 2 * n - 3; m++)
1649  {
1650  double u = 0;
1651  for (k = m + 1 - n; k <= floor((m - 1) / 2.0); k++)
1652  {
1653  int l = m - k;
1654  j = n - 1 - l;
1655  u = u - (a[k + n + 1] - a[l]) * t[j + 1] - b[k + n + 1] * s[j + 1] +
1656  b[l] * s[j + 2];
1657  s[j + 1] = u;
1658  }
1659 
1660  if (m % 2 == 0)
1661  {
1662  k = m / 2;
1663  a[k + n + 1] =
1664  a[k] + (s[j + 1] - b[k + n + 1] * s[j + 2]) / t[j + 2];
1665  }
1666  else
1667  {
1668  k = (m + 1) / 2;
1669  b[k + n + 1] = s[j + 1] / s[j + 2];
1670  }
1671 
1672  // Swap the contents of s and t
1673  double *hold = s;
1674  s = t;
1675  t = hold;
1676  }
1677 
1678  a[2 * n] = a[n - 1] - b[2 * n] * s[1] / t[1];
1679 }

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

◆ laginterp()

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

Definition at line 53 of file Polylib.cpp.

54 {
55  double temp = 1.0;
56  for (int i = 0; i < np; i++)
57  {
58  if (j != i)
59  {
60  temp *= optdiff(z, zj[i]) / (zj[j] - zj[i]);
61  }
62  }
63  return temp;
64 }
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

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

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 
32  m_expn = static_cast<int>(floor(log10(fabs(xl - xr))));
33  m_xln = xl * powl(10.0L, -m_expn) -
34  floor(xl * powl(10.0L,
35  -m_expn)); // substract the digits overlap part
36  m_xrn =
37  xr * powl(10.0L, -m_expn) -
38  floor(xl *
39  powl(10.0L,
40  -m_expn)); // substract the common digits overlap part
41  m_xln =
42  round(m_xln * powl(10.0L, m_digits + m_expn)); // git rid of rubbish
43  m_xrn = round(m_xrn * powl(10.0L, m_digits + m_expn));
44 
45  return powl(10.0L, -m_digits) * (m_xln - m_xrn);
46  }
47  else
48  {
49  return (xl - xr);
50  }
51 }

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

1440 {
1441 
1442  int i;
1443  double apb, apbi, a2b2;
1444 
1445  if (!n)
1446  return;
1447 
1448  // generate normalised terms
1449  apb = alpha + beta;
1450  apbi = 2.0 + apb;
1451 
1452  b[0] = pow(2.0, apb + 1.0) * gammaF(alpha + 1.0) * gammaF(beta + 1.0) /
1453  gammaF(apbi); // MuZero
1454  a[0] = (beta - alpha) / apbi;
1455  b[1] = (4.0 * (1.0 + alpha) * (1.0 + beta) / ((apbi + 1.0) * apbi * apbi));
1456 
1457  a2b2 = beta * beta - alpha * alpha;
1458 
1459  for (i = 1; i < n - 1; i++)
1460  {
1461  apbi = 2.0 * (i + 1) + apb;
1462  a[i] = a2b2 / ((apbi - 2.0) * apbi);
1463  b[i + 1] = (4.0 * (i + 1) * (i + 1 + alpha) * (i + 1 + beta) *
1464  (i + 1 + apb) / ((apbi * apbi - 1) * apbi * apbi));
1465  }
1466 
1467  apbi = 2.0 * n + apb;
1468  a[n - 1] = a2b2 / ((apbi - 2.0) * apbi);
1469 }

References Nektar::LibUtilities::beta, and 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 1497 of file Polylib.cpp.

1498 {
1499  int m, l, iter, i, k;
1500  double s, r, p, g, f, dd, c, b;
1501 
1502  double MuZero = e[0];
1503 
1504  // Renumber the elements of e
1505  for (i = 0; i < n - 1; i++)
1506  {
1507  e[i] = sqrt(e[i + 1]);
1508  }
1509  e[n - 1] = 0.0;
1510 
1511  for (l = 0; l < n; l++)
1512  {
1513  iter = 0;
1514  do
1515  {
1516  for (m = l; m < n - 1; m++)
1517  {
1518  dd = fabs(d[m]) + fabs(d[m + 1]);
1519  if (fabs(e[m]) + dd == dd)
1520  break;
1521  }
1522  if (m != l)
1523  {
1524  if (iter++ == STOP)
1525  {
1526  fprintf(stderr, "triQL: Too many iterations in TQLI");
1527  exit(1);
1528  }
1529  g = (d[l + 1] - d[l]) / (2.0 * e[l]);
1530  r = sqrt((g * g) + 1.0);
1531  g = d[m] - d[l] + e[l] / (g + sign(r, g));
1532  s = c = 1.0;
1533  p = 0.0;
1534  for (i = m - 1; i >= l; i--)
1535  {
1536  f = s * e[i];
1537  b = c * e[i];
1538  if (fabs(f) >= fabs(g))
1539  {
1540  c = g / f;
1541  r = sqrt((c * c) + 1.0);
1542  e[i + 1] = f * r;
1543  c *= (s = 1.0 / r);
1544  }
1545  else
1546  {
1547  s = f / g;
1548  r = sqrt((s * s) + 1.0);
1549  e[i + 1] = g * r;
1550  s *= (c = 1.0 / r);
1551  }
1552  g = d[i + 1] - p;
1553  r = (d[i] - g) * s + 2.0 * c * b;
1554  p = s * r;
1555  d[i + 1] = g + p;
1556  g = c * r - b;
1557 
1558  // Calculate the eigenvectors
1559  for (k = 0; k < n; k++)
1560  {
1561  f = z[k][i + 1];
1562  z[k][i + 1] = s * z[k][i] + c * f;
1563  z[k][i] = c * z[k][i] - s * f;
1564  }
1565  }
1566  d[l] = d[l] - p;
1567  e[l] = g;
1568  e[m] = 0.0;
1569  }
1570  } while (m != l);
1571  }
1572 
1573  // order eigenvalues
1574  // Since we only need the first component of the eigenvectors
1575  // to calcualte the weight, we only swap the first components
1576  for (i = 0; i < n - 1; ++i)
1577  {
1578  k = i;
1579  p = d[i];
1580  for (l = i + 1; l < n; ++l)
1581  if (d[l] < p)
1582  {
1583  k = l;
1584  p = d[l];
1585  }
1586  d[k] = d[i];
1587  d[i] = p;
1588 
1589  double temp = z[0][k];
1590  z[0][k] = z[0][i];
1591  z[0][i] = temp;
1592  }
1593 
1594  // Calculate the weights
1595  for (i = 0; i < n; i++)
1596  {
1597  e[i] = MuZero * z[0][i] * z[0][i];
1598  }
1599 }
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 99 of file Polylib.cpp.

101 {
102  int i;
103  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
104 
105  jacobz(np, z, alpha, beta);
106  jacobd(np, z, w, np, alpha, beta);
107 
108  fac = pow(two, apb + one) * gammaFracGammaF(np + 1, alpha, np + 1, 0.0) *
109  gammaFracGammaF(np + 1, beta, np + 1, apb);
110 
111  for (i = 0; i < np; ++i)
112  w[i] = fac / (w[i] * w[i] * (one - z[i] * z[i]));
113 
114  return;
115 }
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:70

References Nektar::LibUtilities::beta, 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 256 of file Polylib.cpp.

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

References Nektar::LibUtilities::beta, 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 207 of file Polylib.cpp.

209 {
210 
211  if (np == 1)
212  {
213  z[0] = 0.0;
214  w[0] = 2.0;
215  }
216  else if (np == 2)
217  {
218  z[0] = -1.0;
219  z[1] = 1.0;
220 
221  w[0] = 1.0;
222  w[1] = 1.0;
223  }
224  else
225  {
226  int i;
227  double fac, one = 1.0, apb = alpha + beta, two = 2.0;
228 
229  z[0] = -one;
230  z[np - 1] = one;
231  jacobz(np - 2, z + 1, alpha + one, beta + one);
232  jacobfd(np, z, w, NULL, np - 1, alpha, beta);
233 
234  fac = pow(two, apb + 1) * gammaFracGammaF(np, alpha, np, 0.0) *
235  gammaFracGammaF(np, beta, np + 1, apb);
236  fac /= (np - 1);
237 
238  for (i = 0; i < np; ++i)
239  w[i] = fac / (w[i] * w[i]);
240  w[0] *= (beta + one);
241  w[np - 1] *= (alpha + one);
242  }
243 
244  return;
245 }

References Nektar::LibUtilities::beta, 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 127 of file Polylib.cpp.

129 {
130 
131  if (np == 1)
132  {
133  z[0] = 0.0;
134  w[0] = 2.0;
135  }
136  else
137  {
138  int i;
139  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
140 
141  z[0] = -one;
142  jacobz(np - 1, z + 1, alpha, beta + 1);
143  jacobfd(np, z, w, NULL, np - 1, alpha, beta);
144 
145  fac = pow(two, apb) * gammaFracGammaF(np, alpha, np, 0.0) *
146  gammaFracGammaF(np, beta, np + 1, apb);
147  fac /= (beta + np);
148 
149  for (i = 0; i < np; ++i)
150  w[i] = fac * (1 - z[i]) / (w[i] * w[i]);
151  w[0] *= (beta + one);
152  }
153 
154  return;
155 }

References Nektar::LibUtilities::beta, 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 168 of file Polylib.cpp.

170 {
171 
172  if (np == 1)
173  {
174  z[0] = 0.0;
175  w[0] = 2.0;
176  }
177  else
178  {
179  int i;
180  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
181 
182  jacobz(np - 1, z, alpha + 1, beta);
183  z[np - 1] = one;
184  jacobfd(np, z, w, NULL, np - 1, alpha, beta);
185 
186  fac = pow(two, apb) * gammaFracGammaF(np, alpha, np, 0.0) *
187  gammaFracGammaF(np, beta, np + 1, apb);
188  fac /= (alpha + np);
189 
190  for (i = 0; i < np; ++i)
191  w[i] = fac * (1 + z[i]) / (w[i] * w[i]);
192  w[np - 1] *= (alpha + one);
193  }
194 
195  return;
196 }

References Nektar::LibUtilities::beta, 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 434 of file Polylib.cpp.

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

References Nektar::LibUtilities::beta, 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 328 of file Polylib.cpp.

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

References Nektar::LibUtilities::beta, chri1(), JKMatrix(), RecCoeff(), and TriQL().

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

Generated on Tue Aug 23 2022 07:27:24 for Nektar++ by doxygen 1.9.1