Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Functions
Polylib Namespace Reference

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

Functions

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

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

Referenced by zwlk(), and zwrk().

2935  {
2936 
2937 
2938 
2939  double q = ceil(3.0*n/2);
2940 
2941  int size = (int)q+1;
2942 
2943  if(size < n+1)
2944 
2945  {
2946 
2947  fprintf(stderr,"input arrays a and b are too short\n");
2948 
2949  }
2950 
2951  double* r = new double[n+1];
2952 
2953  r[0] = z - a[0];
2954 
2955  r[1] = z - a[1] - b[1]/r[0];
2956 
2957  a0[0] = a[1] + r[1] - r[0];
2958 
2959  b0[0] = -r[0]*b[0];
2960 
2961 
2962 
2963  if(n == 1)
2964 
2965  {
2966 
2967  return;
2968 
2969  }
2970 
2971  int k = 0;
2972 
2973  for(k = 1; k < n; k++)
2974 
2975  {
2976 
2977  r[k+1] = z - a[k+1] - b[k+1]/r[k];
2978 
2979  a0[k] = a[k+1] + r[k+1] - r[k];
2980 
2981  b0[k] = b[k] * r[k]/r[k-1];
2982 
2983  }
2984 
2985 
2986  }
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 937 of file Polylib.cpp.

References jacobd().

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

941  {
942 
943 
944 
945  double one = 1.0, two = 2.0;
946 
947 
948 
949  if (np <= 0){
950 
951  D[0] = 0.0;
952 
953  }
954 
955  else{
956 
957  register int i,j;
958 
959  double *pd;
960 
961 
962 
963  pd = (double *)malloc(np*sizeof(double));
964 
965  jacobd(np,z,pd,np,alpha,beta);
966 
967 
968 
969  for (i = 0; i < np; i++){
970 
971  for (j = 0; j < np; j++){
972 
973 
974 
975  if (i != j)
976 
977  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
978 
979  else
980 
981  D[i*np+j] = (alpha - beta + (alpha + beta + two)*z[j])/
982 
983  (two*(one - z[j]*z[j]));
984 
985  }
986 
987  }
988 
989  free(pd);
990 
991  }
992 
993  return;
994 
995  }
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:2151
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 1221 of file Polylib.cpp.

References gammaF(), and jacobd().

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

1225  {
1226 
1227 
1228 
1229  if (np <= 0){
1230 
1231  D[0] = 0.0;
1232 
1233  }
1234 
1235  else{
1236 
1237  register int i, j;
1238 
1239  double one = 1.0, two = 2.0;
1240 
1241  double *pd;
1242 
1243 
1244 
1245  pd = (double *)malloc(np*sizeof(double));
1246 
1247 
1248 
1249  pd[0] = two*pow(-one,np)*gammaF(np + beta);
1250 
1251  pd[0] /= gammaF(np - one)*gammaF(beta + two);
1252 
1253  jacobd(np-2,z+1,pd+1,np-2,alpha+1,beta+1);
1254 
1255  for(i = 1; i < np-1; ++i) pd[i] *= (one-z[i]*z[i]);
1256 
1257  pd[np-1] = -two*gammaF(np + alpha);
1258 
1259  pd[np-1] /= gammaF(np - one)*gammaF(alpha + two);
1260 
1261 
1262 
1263  for (i = 0; i < np; i++)
1264 
1265  for (j = 0; j < np; j++){
1266 
1267  if (i != j)
1268 
1269  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1270 
1271  else {
1272 
1273  if (j == 0)
1274 
1275  D[i*np+j] = (alpha - (np-1)*(np + alpha + beta))/(two*(beta+ two));
1276 
1277  else if (j == np-1)
1278 
1279  D[i*np+j] =-(beta - (np-1)*(np + alpha + beta))/(two*(alpha+ two));
1280 
1281  else
1282 
1283  D[i*np+j] = (alpha - beta + (alpha + beta)*z[j])/
1284 
1285  (two*(one - z[j]*z[j]));
1286 
1287  }
1288 
1289  }
1290 
1291  free(pd);
1292 
1293  }
1294 
1295 
1296 
1297  return;
1298 
1299  }
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:2151
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2205
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 1023 of file Polylib.cpp.

References gammaF(), and jacobd().

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

1027  {
1028 
1029 
1030 
1031  if (np <= 0){
1032 
1033  D[0] = 0.0;
1034 
1035  }
1036 
1037  else{
1038 
1039  register int i, j;
1040 
1041  double one = 1.0, two = 2.0;
1042 
1043  double *pd;
1044 
1045 
1046 
1047  pd = (double *)malloc(np*sizeof(double));
1048 
1049 
1050 
1051  pd[0] = pow(-one,np-1)*gammaF(np+beta+one);
1052 
1053  pd[0] /= gammaF(np)*gammaF(beta+two);
1054 
1055  jacobd(np-1,z+1,pd+1,np-1,alpha,beta+1);
1056 
1057  for(i = 1; i < np; ++i) pd[i] *= (1+z[i]);
1058 
1059 
1060 
1061  for (i = 0; i < np; i++)
1062 
1063  for (j = 0; j < np; j++){
1064 
1065  if (i != j)
1066 
1067  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1068 
1069  else {
1070 
1071  if(j == 0)
1072 
1073  D[i*np+j] = -(np + alpha + beta + one)*(np - one)/
1074 
1075  (two*(beta + two));
1076 
1077  else
1078 
1079  D[i*np+j] = (alpha - beta + one + (alpha + beta + one)*z[j])/
1080 
1081  (two*(one - z[j]*z[j]));
1082 
1083  }
1084 
1085  }
1086 
1087  free(pd);
1088 
1089  }
1090 
1091 
1092 
1093  return;
1094 
1095  }
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:2151
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2205
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 1121 of file Polylib.cpp.

References gammaF(), and jacobd().

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

1125  {
1126 
1127 
1128 
1129  if (np <= 0){
1130 
1131  D[0] = 0.0;
1132 
1133  }
1134 
1135  else{
1136 
1137  register int i, j;
1138 
1139  double one = 1.0, two = 2.0;
1140 
1141  double *pd;
1142 
1143 
1144 
1145  pd = (double *)malloc(np*sizeof(double));
1146 
1147 
1148 
1149 
1150 
1151  jacobd(np-1,z,pd,np-1,alpha+1,beta);
1152 
1153  for(i = 0; i < np-1; ++i) pd[i] *= (1-z[i]);
1154 
1155  pd[np-1] = -gammaF(np+alpha+one);
1156 
1157  pd[np-1] /= gammaF(np)*gammaF(alpha+two);
1158 
1159 
1160 
1161  for (i = 0; i < np; i++)
1162 
1163  for (j = 0; j < np; j++){
1164 
1165  if (i != j)
1166 
1167  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1168 
1169  else {
1170 
1171  if(j == np-1)
1172 
1173  D[i*np+j] = (np + alpha + beta + one)*(np - one)/
1174 
1175  (two*(alpha + two));
1176 
1177  else
1178 
1179  D[i*np+j] = (alpha - beta - one + (alpha + beta + one)*z[j])/
1180 
1181  (two*(one - z[j]*z[j]));
1182 
1183  }
1184 
1185  }
1186 
1187  free(pd);
1188 
1189  }
1190 
1191 
1192 
1193  return;
1194 
1195  }
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:2151
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2205
double Polylib::gammaF ( const double  x)

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

values and halves.

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

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

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

Definition at line 2205 of file Polylib.cpp.

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

2205  {
2206 
2207  double gamma = 1.0;
2208 
2209 
2210 
2211  if (x == -0.5) gamma = -2.0*sqrt(M_PI);
2212 
2213  else if (!x) return gamma;
2214 
2215  else if ((x-(int)x) == 0.5){
2216 
2217  int n = (int) x;
2218 
2219  double tmp = x;
2220 
2221 
2222 
2223  gamma = sqrt(M_PI);
2224 
2225  while(n--){
2226 
2227  tmp -= 1.0;
2228 
2229  gamma *= tmp;
2230 
2231  }
2232 
2233  }
2234 
2235  else if ((x-(int)x) == 0.0){
2236 
2237  int n = (int) x;
2238 
2239  double tmp = x;
2240 
2241 
2242 
2243  while(--n){
2244 
2245  tmp -= 1.0;
2246 
2247  gamma *= tmp;
2248 
2249  }
2250 
2251  }
2252 
2253  else
2254 
2255  fprintf(stderr,"%lf is not of integer or half order\n",x);
2256 
2257  return gamma;
2258 
2259  }
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 1345 of file Polylib.cpp.

References EPS, jacobd(), jacobfd(), and CellMLToNektar.cellml_metadata::p.

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

1349  {
1350 
1351 
1352 
1353  double zi, dz, p, pd, h;
1354 
1355 
1356 
1357  zi = *(zgj+i);
1358 
1359  dz = z - zi;
1360 
1361  if (fabs(dz) < EPS) return 1.0;
1362 
1363 
1364 
1365  jacobd (1, &zi, &pd , np, alpha, beta);
1366 
1367  jacobfd(1, &z , &p, NULL , np, alpha, beta);
1368 
1369  h = p/(pd*dz);
1370 
1371 
1372 
1373  return h;
1374 
1375  }
CellMLToNektar.cellml_metadata.p
tuple p
Definition: cellml_metadata.py:350
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:2151
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:23
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:1951
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 1589 of file Polylib.cpp.

References EPS, jacobd(), jacobfd(), and CellMLToNektar.cellml_metadata::p.

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

1593  {
1594 
1595  double one = 1., two = 2.;
1596 
1597  double zi, dz, p, pd, h;
1598 
1599 
1600 
1601  zi = *(zglj+i);
1602 
1603  dz = z - zi;
1604 
1605  if (fabs(dz) < EPS) return 1.0;
1606 
1607 
1608 
1609  jacobfd(1, &zi, &p , NULL, np-2, alpha + one, beta + one);
1610 
1611  // need to use this routine in caes z = -1 or 1
1612 
1613  jacobd (1, &zi, &pd, np-2, alpha + one, beta + one);
1614 
1615  h = (one - zi*zi)*pd - two*zi*p;
1616 
1617  jacobfd(1, &z, &p, NULL, np-2, alpha + one, beta + one);
1618 
1619  h = (one - z*z)*p/(h*dz);
1620 
1621 
1622 
1623  return h;
1624 
1625  }
CellMLToNektar.cellml_metadata.p
tuple p
Definition: cellml_metadata.py:350
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:2151
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:23
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:1951
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 1421 of file Polylib.cpp.

References EPS, jacobd(), jacobfd(), and CellMLToNektar.cellml_metadata::p.

Referenced by Imgrjm().

1425  {
1426 
1427 
1428 
1429  double zi, dz, p, pd, h;
1430 
1431 
1432 
1433  zi = *(zgrj+i);
1434 
1435  dz = z - zi;
1436 
1437  if (fabs(dz) < EPS) return 1.0;
1438 
1439 
1440 
1441  jacobfd (1, &zi, &p , NULL, np-1, alpha, beta + 1);
1442 
1443  // need to use this routine in caes zi = -1 or 1
1444 
1445  jacobd (1, &zi, &pd, np-1, alpha, beta + 1);
1446 
1447  h = (1.0 + zi)*pd + p;
1448 
1449  jacobfd (1, &z, &p, NULL, np-1, alpha, beta + 1);
1450 
1451  h = (1.0 + z )*p/(h*dz);
1452 
1453 
1454 
1455  return h;
1456 
1457  }
CellMLToNektar.cellml_metadata.p
tuple p
Definition: cellml_metadata.py:350
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:2151
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:23
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:1951
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 1505 of file Polylib.cpp.

References EPS, jacobd(), jacobfd(), and CellMLToNektar.cellml_metadata::p.

Referenced by Imgrjp().

1509  {
1510 
1511 
1512 
1513  double zi, dz, p, pd, h;
1514 
1515 
1516 
1517  zi = *(zgrj+i);
1518 
1519  dz = z - zi;
1520 
1521  if (fabs(dz) < EPS) return 1.0;
1522 
1523 
1524 
1525  jacobfd (1, &zi, &p , NULL, np-1, alpha+1, beta );
1526 
1527  // need to use this routine in caes z = -1 or 1
1528 
1529  jacobd (1, &zi, &pd, np-1, alpha+1, beta );
1530 
1531  h = (1.0 - zi)*pd - p;
1532 
1533  jacobfd (1, &z, &p, NULL, np-1, alpha+1, beta);
1534 
1535  h = (1.0 - z )*p/(h*dz);
1536 
1537 
1538 
1539  return h;
1540 
1541  }
CellMLToNektar.cellml_metadata.p
tuple p
Definition: cellml_metadata.py:350
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:2151
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:23
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:1951
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 2999 of file Polylib.cpp.

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

3000  {
3001  std::complex<Nektar::NekDouble> z (1.0,0.0);
3002  std::complex<Nektar::NekDouble> zbes (1.0,0.0);
3003  std::complex<Nektar::NekDouble> zarg;
3004  Nektar::NekDouble tol = 1e-15;
3005  int maxit = 10000;
3006  int i = 1;
3007 
3008  zarg = -0.25*y*y;
3009 
3010  while (abs(z) > tol && i <= maxit){
3011  z = z*(1.0/i/(i+n)*zarg);
3012  if (abs(z) <= tol) break;
3013  zbes = zbes + z;
3014  i++;
3015  }
3016  zarg = 0.5*y;
3017  for (i=1;i<=n;i++){
3018  zbes = zbes*zarg;
3019  }
3020  return zbes;
3021 
3022  }
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
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 1657 of file Polylib.cpp.

References hgj().

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

1659  {
1660 
1661  double zp;
1662 
1663  register int i, j;
1664 
1665 
1666 
1667  for (i = 0; i < nz; ++i) {
1668 
1669  for (j = 0; j < mz; ++j)
1670 
1671  {
1672 
1673  zp = zm[j];
1674 
1675  im [i*mz+j] = hgj(i, zp, zgj, nz, alpha, beta);
1676 
1677  }
1678 
1679  }
1680 
1681 
1682 
1683  return;
1684 
1685  }
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.
Definition: Polylib.cpp:1345
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 1837 of file Polylib.cpp.

References hglj().

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

1841  {
1842 
1843  double zp;
1844 
1845  register int i, j;
1846 
1847 
1848 
1849  for (i = 0; i < nz; i++) {
1850 
1851  for (j = 0; j < mz; j++)
1852 
1853  {
1854 
1855  zp = zm[j];
1856 
1857  im[i*mz+j] = hglj(i, zp, zglj, nz, alpha, beta);
1858 
1859  }
1860 
1861  }
1862 
1863 
1864 
1865  return;
1866 
1867  }
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.
Definition: Polylib.cpp:1589
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 1715 of file Polylib.cpp.

References hgrjm().

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

1719  {
1720 
1721  double zp;
1722 
1723  register int i, j;
1724 
1725 
1726 
1727  for (i = 0; i < nz; i++) {
1728 
1729  for (j = 0; j < mz; j++)
1730 
1731  {
1732 
1733  zp = zm[j];
1734 
1735  im [i*mz+j] = hgrjm(i, zp, zgrj, nz, alpha, beta);
1736 
1737  }
1738 
1739  }
1740 
1741 
1742 
1743  return;
1744 
1745  }
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.
Definition: Polylib.cpp:1421
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 1775 of file Polylib.cpp.

References hgrjp().

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

1779  {
1780 
1781  double zp;
1782 
1783  register int i, j;
1784 
1785 
1786 
1787  for (i = 0; i < nz; i++) {
1788 
1789  for (j = 0; j < mz; j++)
1790 
1791  {
1792 
1793  zp = zm[j];
1794 
1795  im [i*mz+j] = hgrjp(i, zp, zgrj, nz, alpha, beta);
1796 
1797  }
1798 
1799  }
1800 
1801 
1802 
1803  return;
1804 
1805  }
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.
Definition: Polylib.cpp:1505
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 2151 of file Polylib.cpp.

References jacobfd().

Referenced by Dgj(), Dglj(), Dgrjm(), Dgrjp(), hgj(), hglj(), hgrjm(), hgrjp(), 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().

2155  {
2156 
2157  register int i;
2158 
2159  double one = 1.0;
2160 
2161  if(n == 0)
2162 
2163  for(i = 0; i < np; ++i) polyd[i] = 0.0;
2164 
2165  else{
2166 
2167  //jacobf(np,z,polyd,n-1,alpha+one,beta+one);
2168 
2169  jacobfd(np,z,polyd,NULL,n-1,alpha+one,beta+one);
2170 
2171  for(i = 0; i < np; ++i) polyd[i] *= 0.5*(alpha + beta + (double)n + one);
2172 
2173  }
2174 
2175  return;
2176 
2177  }
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:1951
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 1951 of file Polylib.cpp.

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

1953  {
1954 
1955  register int i;
1956 
1957  double zero = 0.0, one = 1.0, two = 2.0;
1958 
1959 
1960 
1961  if(!np)
1962 
1963  return;
1964 
1965 
1966 
1967  if(n == 0){
1968 
1969  if(poly_in)
1970 
1971  for(i = 0; i < np; ++i)
1972 
1973  poly_in[i] = one;
1974 
1975  if(polyd)
1976 
1977  for(i = 0; i < np; ++i)
1978 
1979  polyd[i] = zero;
1980 
1981  }
1982 
1983  else if (n == 1){
1984 
1985  if(poly_in)
1986 
1987  for(i = 0; i < np; ++i)
1988 
1989  poly_in[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
1990 
1991  if(polyd)
1992 
1993  for(i = 0; i < np; ++i)
1994 
1995  polyd[i] = 0.5*(alpha + beta + two);
1996 
1997  }
1998 
1999  else{
2000 
2001  register int k;
2002 
2003  double a1,a2,a3,a4;
2004 
2005  double two = 2.0, apb = alpha + beta;
2006 
2007  double *poly, *polyn1,*polyn2;
2008 
2009 
2010 
2011  if(poly_in){ // switch for case of no poynomial function return
2012 
2013  polyn1 = (double *)malloc(2*np*sizeof(double));
2014 
2015  polyn2 = polyn1+np;
2016 
2017  poly = poly_in;
2018 
2019  }
2020 
2021  else{
2022 
2023  polyn1 = (double *)malloc(3*np*sizeof(double));
2024 
2025  polyn2 = polyn1+np;
2026 
2027  poly = polyn2+np;
2028 
2029  }
2030 
2031 
2032 
2033  for(i = 0; i < np; ++i){
2034 
2035  polyn2[i] = one;
2036 
2037  polyn1[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
2038 
2039  }
2040 
2041 
2042 
2043  for(k = 2; k <= n; ++k){
2044 
2045  a1 = two*k*(k + apb)*(two*k + apb - two);
2046 
2047  a2 = (two*k + apb - one)*(alpha*alpha - beta*beta);
2048 
2049  a3 = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb);
2050 
2051  a4 = two*(k + alpha - one)*(k + beta - one)*(two*k + apb);
2052 
2053 
2054 
2055  a2 /= a1;
2056 
2057  a3 /= a1;
2058 
2059  a4 /= a1;
2060 
2061 
2062 
2063  for(i = 0; i < np; ++i){
2064 
2065  poly [i] = (a2 + a3*z[i])*polyn1[i] - a4*polyn2[i];
2066 
2067  polyn2[i] = polyn1[i];
2068 
2069  polyn1[i] = poly [i];
2070 
2071  }
2072 
2073  }
2074 
2075 
2076 
2077  if(polyd){
2078 
2079  a1 = n*(alpha - beta);
2080 
2081  a2 = n*(two*n + alpha + beta);
2082 
2083  a3 = two*(n + alpha)*(n + beta);
2084 
2085  a4 = (two*n + alpha + beta);
2086 
2087  a1 /= a4; a2 /= a4; a3 /= a4;
2088 
2089 
2090 
2091  // note polyn2 points to polyn1 at end of poly iterations
2092 
2093  for(i = 0; i < np; ++i){
2094 
2095  polyd[i] = (a1- a2*z[i])*poly[i] + a3*polyn2[i];
2096 
2097  polyd[i] /= (one - z[i]*z[i]);
2098 
2099  }
2100 
2101  }
2102 
2103 
2104 
2105  free(polyn1);
2106 
2107  }
2108 
2109 
2110 
2111  return;
2112 
2113  }
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 2279 of file Polylib.cpp.

References EPS, jacobfd(), and STOP.

2281  {
2282 
2283  register int i,j,k;
2284 
2285  double dth = M_PI/(2.0*(double)n);
2286 
2287  double poly,pder,rlast=0.0;
2288 
2289  double sum,delr,r;
2290 
2291  double one = 1.0, two = 2.0;
2292 
2293 
2294 
2295  if(!n)
2296 
2297  return;
2298 
2299 
2300 
2301  for(k = 0; k < n; ++k){
2302 
2303  r = -cos((two*(double)k + one) * dth);
2304 
2305  if(k) r = 0.5*(r + rlast);
2306 
2307 
2308 
2309  for(j = 1; j < STOP; ++j){
2310 
2311  jacobfd(1,&r,&poly, &pder, n, alpha, beta);
2312 
2313 
2314 
2315  for(i = 0, sum = 0.0; i < k; ++i) sum += one/(r - z[i]);
2316 
2317 
2318 
2319  delr = -poly / (pder - sum * poly);
2320 
2321  r += delr;
2322 
2323  if( fabs(delr) < EPS ) break;
2324 
2325  }
2326 
2327  z[k] = r;
2328 
2329  rlast = r;
2330 
2331  }
2332 
2333  return;
2334 
2335  }
STOP
#define STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
Definition: Polylib.cpp:19
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:23
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:1951
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 2389 of file Polylib.cpp.

References RecCoeff(), and TriQL().

2391  {
2392 
2393 
2394 
2395  int i,j;
2396 
2397  RecCoeff(n,a,b,alpha,beta);
2398 
2399 
2400 
2401  double **z = new double*[n];
2402 
2403  for(i = 0; i < n; i++)
2404 
2405  {
2406 
2407  z[i] = new double[n];
2408 
2409  for(j = 0; j < n; j++)
2410 
2411  {
2412 
2413  z[i][j] = 0.0;
2414 
2415  }
2416 
2417  }
2418 
2419  for(i = 0; i < n; i++)
2420 
2421  {
2422 
2423  z[i][i] = 1.0;
2424 
2425  }
2426 
2427 
2428 
2429  // find eigenvalues and eigenvectors
2430 
2431  TriQL(n, a, b,z);
2432 
2433 
2434 
2435  return;
2436 
2437  }
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
Definition: Polylib.cpp:2449
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:2565
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 2769 of file Polylib.cpp.

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

2771  {
2772 
2773  int i,j,k,m;
2774 
2775  // Working storage
2776 
2777  int size = (int)floor(n/2.0)+2;
2778 
2779  double *s = new double[size];
2780 
2781  double *t = new double[size];
2782 
2783 
2784 
2785  // Initialize s and t to zero
2786 
2787  for(i = 0; i < size; i++)
2788 
2789  {
2790 
2791  s[i] = 0.0;
2792 
2793  t[i] = 0.0;
2794 
2795  }
2796 
2797 
2798 
2799  t[1] = b[n+1];
2800 
2801  for(m = 0; m <= n-2; m++)
2802 
2803  {
2804 
2805  double u = 0.0;
2806 
2807  for(k = (int)floor((m+1)/2.0); k >= 0; k--)
2808 
2809  {
2810 
2811  int l = m-k;
2812 
2813  u = u+(a[k+n+1]-a[l])*t[k+1] + b[k+n+1]*s[k] - b[l]*s[k+1];
2814 
2815  s[k+1] = u;
2816 
2817  }
2818 
2819 
2820 
2821  // Swap the contents of s and t
2822 
2823  double *hold = s;
2824 
2825  s = t;
2826 
2827  t = hold;
2828 
2829  }
2830 
2831 
2832 
2833 
2834 
2835  for(j = (int)floor(n/2.0); j >= 0; j--)
2836 
2837  {
2838 
2839  s[j+1] = s[j];
2840 
2841  }
2842 
2843 
2844 
2845  for(m = n-1; m <= 2*n-3; m++)
2846 
2847  {
2848 
2849  double u = 0;
2850 
2851  for(k = m+1-n; k <= floor((m-1)/2.0); k++)
2852 
2853  {
2854 
2855  int l = m-k;
2856 
2857  j = n-1-l;
2858 
2859  u = u-(a[k+n+1]-a[l])*t[j+1] - b[k+n+1]*s[j+1] + b[l]*s[j+2];
2860 
2861  s[j+1] = u;
2862 
2863  }
2864 
2865 
2866 
2867  if(m%2 == 0)
2868 
2869  {
2870 
2871  k = m/2;
2872 
2873  a[k+n+1] = a[k] + (s[j+1]-b[k+n+1]*s[j+2])/t[j+2];
2874 
2875 
2876 
2877  }else
2878 
2879  {
2880 
2881  k = (m+1)/2;
2882 
2883  b[k+n+1] = s[j+1]/s[j+2];
2884 
2885  }
2886 
2887 
2888 
2889 
2890 
2891  // Swap the contents of s and t
2892 
2893  double *hold = s;
2894 
2895  s = t;
2896 
2897  t = hold;
2898 
2899  }
2900 
2901 
2902 
2903  a[2*n ] = a[n-1]-b[2*n]*s[1]/t[1];
2904 
2905 
2906 
2907  }
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 2449 of file Polylib.cpp.

References gammaF().

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

2451  {
2452 
2453 
2454 
2455  int i;
2456 
2457  double apb, apbi,a2b2;
2458 
2459 
2460 
2461  if(!n)
2462 
2463  return;
2464 
2465 
2466 
2467  // generate normalised terms
2468 
2469  apb = alpha + beta;
2470 
2471  apbi = 2.0 + apb;
2472 
2473 
2474 
2475  b[0] = pow(2.0,apb+1.0)*gammaF(alpha+1.0)*gammaF(beta+1.0)/gammaF(apbi); //MuZero
2476 
2477  a[0] = (beta-alpha)/apbi;
2478 
2479  b[1] = (4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));
2480 
2481 
2482 
2483  a2b2 = beta*beta-alpha*alpha;
2484 
2485 
2486 
2487  for(i = 1; i < n-1; i++){
2488 
2489  apbi = 2.0*(i+1) + apb;
2490 
2491  a[i] = a2b2/((apbi-2.0)*apbi);
2492 
2493  b[i+1] = (4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
2494 
2495  ((apbi*apbi-1)*apbi*apbi));
2496 
2497  }
2498 
2499 
2500 
2501  apbi = 2.0*n + apb;
2502 
2503  a[n-1] = a2b2/((apbi-2.0)*apbi);
2504 
2505 
2506 
2507  }
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2205
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 2565 of file Polylib.cpp.

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

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

2565  {
2566 
2567  int m,l,iter,i,k;
2568 
2569  double s,r,p,g,f,dd,c,b;
2570 
2571 
2572 
2573  double MuZero = e[0];
2574 
2575 
2576 
2577  // Renumber the elements of e
2578 
2579  for(i = 0; i < n-1; i++)
2580 
2581  {
2582 
2583  e[i] = sqrt(e[i+1]);
2584 
2585  }
2586 
2587  e[n-1] = 0.0;
2588 
2589 
2590 
2591 
2592 
2593  for (l=0;l<n;l++) {
2594 
2595  iter=0;
2596 
2597  do {
2598 
2599  for (m=l;m<n-1;m++) {
2600 
2601  dd=fabs(d[m])+fabs(d[m+1]);
2602 
2603  if (fabs(e[m])+dd == dd) break;
2604 
2605  }
2606 
2607  if (m != l) {
2608 
2609  if (iter++ == STOP){
2610 
2611  fprintf(stderr,"triQL: Too many iterations in TQLI");
2612 
2613  exit(1);
2614 
2615  }
2616 
2617  g=(d[l+1]-d[l])/(2.0*e[l]);
2618 
2619  r=sqrt((g*g)+1.0);
2620 
2621  g=d[m]-d[l]+e[l]/(g+sign(r,g));
2622 
2623  s=c=1.0;
2624 
2625  p=0.0;
2626 
2627  for (i=m-1;i>=l;i--) {
2628 
2629  f=s*e[i];
2630 
2631  b=c*e[i];
2632 
2633  if (fabs(f) >= fabs(g)) {
2634 
2635  c=g/f;
2636 
2637  r=sqrt((c*c)+1.0);
2638 
2639  e[i+1]=f*r;
2640 
2641  c *= (s=1.0/r);
2642 
2643  } else {
2644 
2645  s=f/g;
2646 
2647  r=sqrt((s*s)+1.0);
2648 
2649  e[i+1]=g*r;
2650 
2651  s *= (c=1.0/r);
2652 
2653  }
2654 
2655  g=d[i+1]-p;
2656 
2657  r=(d[i]-g)*s+2.0*c*b;
2658 
2659  p=s*r;
2660 
2661  d[i+1]=g+p;
2662 
2663  g=c*r-b;
2664 
2665 
2666 
2667  // Calculate the eigenvectors
2668 
2669  for(k = 0; k < n; k++)
2670 
2671  {
2672 
2673  f = z[k][i+1];
2674 
2675  z[k][i+1] = s*z[k][i] + c*f;
2676 
2677  z[k][i] = c*z[k][i] - s*f;
2678 
2679  }
2680 
2681 
2682 
2683  }
2684 
2685  d[l]=d[l]-p;
2686 
2687  e[l]=g;
2688 
2689  e[m]=0.0;
2690 
2691  }
2692 
2693  } while (m != l);
2694 
2695  }
2696 
2697 
2698 
2699  // order eigenvalues
2700 
2701  // Since we only need the first component of the eigenvectors
2702 
2703  // to calcualte the weight, we only swap the first components
2704 
2705  for(i = 0; i < n-1; ++i){
2706 
2707  k = i;
2708 
2709  p = d[i];
2710 
2711  for(l = i+1; l < n; ++l)
2712 
2713  if (d[l] < p) {
2714 
2715  k = l;
2716 
2717  p = d[l];
2718 
2719  }
2720 
2721  d[k] = d[i];
2722 
2723  d[i] = p;
2724 
2725 
2726 
2727  double temp = z[0][k];
2728 
2729  z[0][k] = z[0][i];
2730 
2731  z[0][i] = temp;
2732 
2733  }
2734 
2735 
2736 
2737  // Calculate the weights
2738 
2739  for(i =0 ; i < n; i++)
2740 
2741  {
2742 
2743  e[i] = MuZero*z[0][i]*z[0][i];
2744 
2745  }
2746 
2747  }
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:27
CellMLToNektar.cellml_metadata.p
tuple p
Definition: cellml_metadata.py:350
STOP
#define STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
Definition: Polylib.cpp:19
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 107 of file Polylib.cpp.

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

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

111  {
112 
113  register int i;
114 
115  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
116 
117 
118 
119  jacobz (np,z,alpha,beta);
120 
121  jacobd (np,z,w,np,alpha,beta);
122 
123 
124 
125  fac = pow(two,apb + one)*gammaF(alpha + np + one)*gammaF(beta + np + one);
126 
127  fac /= gammaF(np + one)*gammaF(apb + np + one);
128 
129 
130 
131  for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]*(one-z[i]*z[i]));
132 
133 
134 
135  return;
136 
137  }
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:2151
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:45
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2205
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 385 of file Polylib.cpp.

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

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

389  {
390 
391 
392 
393  int np = (npt-1)/2;
394 
395 
396 
397  int i,j;
398 
399 
400 
401  // number of kronrod points associated with the np gauss rule
402 
403  int kpoints = 2*np + 1;
404 
405 
406 
407  // Define the number of required recurrence coefficents
408 
409  int ncoeffs = (int)floor(3.0*(np+1)/2);
410 
411 
412 
413  // Define arrays for the recurrence coefficients
414 
415  // We will use these arrays for the Kronrod results too, hence the
416 
417  // reason for the size of the arrays
418 
419  double *a = new double[kpoints];
420 
421  double *b = new double[kpoints];
422 
423 
424 
425  // Initialize a and b to zero
426 
427  for(i = 0; i < kpoints; i++)
428 
429  {
430 
431  a[i] = 0.0;
432 
433  b[i] = 0.0;
434 
435  }
436 
437 
438 
439  // Call the routine to calculate the recurrence coefficients
440 
441  RecCoeff(ncoeffs,a,b,alpha,beta);
442 
443 
444 
445  // Call the routine to caluclate the jacobi-Kronrod matrix
446 
447  JKMatrix(np,a,b);
448 
449 
450 
451  // Set up the identity matrix
452 
453  double** zmatrix = new double*[kpoints];
454 
455  for(i = 0; i < kpoints; i++)
456 
457  {
458 
459  zmatrix[i] = new double[kpoints];
460 
461  for(j = 0; j < kpoints; j++)
462 
463  {
464 
465  zmatrix[i][j] = 0.0;
466 
467  }
468 
469  }
470 
471  for(i = 0; i < kpoints; i++)
472 
473  {
474 
475  zmatrix[i][i] = 1.0;
476 
477  }
478 
479 
480 
481  // Calculte the points and weights
482 
483  TriQL(kpoints, a, b, zmatrix);
484 
485 
486 
487  for(i = 0; i < kpoints; i++)
488 
489  {
490 
491  z[i] = a[i];
492 
493  w[i] = b[i];
494 
495  }
496 
497 
498 
499 
500 
501  }
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
Definition: Polylib.cpp:2449
Polylib::JKMatrix
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the.
Definition: Polylib.cpp:2769
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:2565
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 311 of file Polylib.cpp.

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

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

315  {
316 
317 
318 
319  if( np == 1 ){
320 
321  z[0] = 0.0;
322 
323  w[0] = 2.0;
324 
325  }
326 
327  else{
328 
329  register int i;
330 
331  double fac, one = 1.0, apb = alpha + beta, two = 2.0;
332 
333 
334 
335  z[0] = -one;
336 
337  z[np-1] = one;
338 
339  jacobz (np-2,z + 1,alpha + one,beta + one);
340 
341  jacobfd (np,z,w,NULL,np-1,alpha,beta);
342 
343 
344 
345  fac = pow(two,apb + 1)*gammaF(alpha + np)*gammaF(beta + np);
346 
347  fac /= (np-1)*gammaF(np)*gammaF(alpha + beta + np + one);
348 
349 
350 
351  for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]);
352 
353  w[0] *= (beta + one);
354 
355  w[np-1] *= (alpha + one);
356 
357  }
358 
359 
360 
361  return;
362 
363  }
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:45
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:1951
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2205
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 163 of file Polylib.cpp.

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

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

167  {
168 
169 
170 
171  if(np == 1){
172 
173  z[0] = 0.0;
174 
175  w[0] = 2.0;
176 
177  }
178 
179  else{
180 
181  register int i;
182 
183  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
184 
185 
186 
187  z[0] = -one;
188 
189  jacobz (np-1,z+1,alpha,beta+1);
190 
191  jacobfd (np,z,w,NULL,np-1,alpha,beta);
192 
193 
194 
195  fac = pow(two,apb)*gammaF(alpha + np)*gammaF(beta + np);
196 
197  fac /= gammaF(np)*(beta + np)*gammaF(apb + np + 1);
198 
199 
200 
201  for(i = 0; i < np; ++i) w[i] = fac*(1-z[i])/(w[i]*w[i]);
202 
203  w[0] *= (beta + one);
204 
205  }
206 
207 
208 
209  return;
210 
211  }
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:45
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:1951
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2205
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 239 of file Polylib.cpp.

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

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

243  {
244 
245 
246 
247  if(np == 1){
248 
249  z[0] = 0.0;
250 
251  w[0] = 2.0;
252 
253  }
254 
255  else{
256 
257  register int i;
258 
259  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
260 
261 
262 
263  jacobz (np-1,z,alpha+1,beta);
264 
265  z[np-1] = one;
266 
267  jacobfd (np,z,w,NULL,np-1,alpha,beta);
268 
269 
270 
271  fac = pow(two,apb)*gammaF(alpha + np)*gammaF(beta + np);
272 
273  fac /= gammaF(np)*(alpha + np)*gammaF(apb + np + 1);
274 
275 
276 
277  for(i = 0; i < np; ++i) w[i] = fac*(1+z[i])/(w[i]*w[i]);
278 
279  w[np-1] *= (alpha + one);
280 
281  }
282 
283 
284 
285  return;
286 
287  }
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:45
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:1951
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2205
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 711 of file Polylib.cpp.

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

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

715  {
716 
717 
718 
719  int np = (npt+1)/2;
720 
721 
722 
723  if(np < 4)
724 
725  {
726 
727  fprintf (stderr,"too few points in formula\n");
728 
729  return;
730 
731  }
732 
733 
734 
735  double endl = -1;
736 
737  double endr = 1;
738 
739  int i,j;
740 
741 
742 
743  // number of kronrod points associated with the np gauss rule
744 
745  int kpoints = 2*np-1;
746 
747 
748 
749  // Define the number of required recurrence coefficents
750 
751  int ncoeffs = (int)ceil(3.0*np/2)-1;
752 
753 
754 
755  // Define arrays for the recurrence coefficients
756 
757  double *a = new double[ncoeffs+1];
758 
759  double *b = new double[ncoeffs+1];
760 
761 
762 
763  // Initialize a and b to zero
764 
765  for(i = 0; i < ncoeffs+1; i++)
766 
767  {
768 
769  a[i] = 0.0;
770 
771  b[i] = 0.0;
772 
773  }
774 
775 
776 
777  // Call the routine to calculate the recurrence coefficients
778 
779  RecCoeff(ncoeffs,a,b,alpha,beta);
780 
781 
782 
783 
784 
785  double* a0 = new double[ncoeffs];
786 
787  double* b0 = new double[ncoeffs];
788 
789 
790 
791  chri1(ncoeffs,a,b,a0,b0,endl);
792 
793 
794 
795  double* a1 = new double[ncoeffs-1];
796 
797  double* b1 = new double[ncoeffs-1];
798 
799 
800 
801  chri1(ncoeffs-1,a0,b0,a1,b1,endr);
802 
803 
804 
805 
806 
807  double s = b1[0]/fabs(b1[0]);
808 
809  b1[0] = s*b1[0];
810 
811 
812 
813  // Finding the 2*np-1 gauss-kronrod points
814 
815  double* z1 = new double[2*np-3];
816 
817  double* w1 = new double[2*np-3];
818 
819  for(i = 0; i < ncoeffs; i++)
820 
821  {
822 
823  z1[i] = a1[i];
824 
825  w1[i] = b1[i];
826 
827  }
828 
829  JKMatrix(np-2,z1,w1);
830 
831  // Set up the identity matrix
832 
833  double** zmatrix = new double*[2*np-3];
834 
835  for(i = 0; i < 2*np-3; i++)
836 
837  {
838 
839  zmatrix[i] = new double[2*np-3];
840 
841  for(j = 0; j < 2*np-3; j++)
842 
843  {
844 
845  zmatrix[i][j] = 0.0;
846 
847  }
848 
849  }
850 
851  for(i = 0; i < 2*np-3; i++)
852 
853  {
854 
855  zmatrix[i][i] = 1.0;
856 
857  }
858 
859 
860 
861  // Calculate the points and weights
862 
863  TriQL(2*np-3, z1, w1, zmatrix);
864 
865 
866 
867  double sumW1 = 0.0;
868 
869  double sumW1Z1 = 0.0;
870 
871  for(i = 0; i < 2*np-3; i++)
872 
873  {
874 
875  w1[i] = s*w1[i]/(z1[i]-endl)/(z1[i]-endr);
876 
877  sumW1 += w1[i];
878 
879  sumW1Z1 += z1[i]*w1[i];
880 
881  }
882 
883 
884 
885  double c0 = b[0]-sumW1;
886 
887  double c1 = a[0]*b[0]-sumW1Z1;
888 
889 
890 
891  z[0] = endl;
892 
893  z[2*np-2] = endr;
894 
895  w[0] = (c0*endr-c1)/(endr-endl);
896 
897  w[2*np-2] = (c1-c0*endl)/(endr-endl);
898 
899 
900 
901  for(i = 1; i < kpoints-1; i++)
902 
903  {
904 
905  z[i] = z1[i-1];
906 
907  w[i] = w1[i-1];
908 
909  }
910 
911  }
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
Definition: Polylib.cpp:2449
Polylib::chri1
void chri1(int, double *, double *, double *, double *, double)
Given a weight function through the first n+1.
Definition: Polylib.cpp:2931
Polylib::JKMatrix
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the.
Definition: Polylib.cpp:2769
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:2565
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 519 of file Polylib.cpp.

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

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

523  {
524 
525 
526 
527  int np = npt/2;
528 
529 
530 
531  if(np < 2)
532 
533  {
534 
535  fprintf(stderr,"too few points in formula\n");
536 
537  return;
538 
539  }
540 
541 
542 
543  double end0 = -1;
544 
545 
546 
547  int i,j;
548 
549 
550 
551  // number of kronrod points associated with the np gauss rule
552 
553  int kpoints = 2*np;
554 
555 
556 
557  // Define the number of required recurrence coefficents
558 
559  int ncoeffs = (int)ceil(3.0*np/2);
560 
561 
562 
563  // Define arrays for the recurrence coefficients
564 
565  double *a = new double[ncoeffs+1];
566 
567  double *b = new double[ncoeffs+1];
568 
569 
570 
571  // Initialize a and b to zero
572 
573  for(i = 0; i < ncoeffs+1; i++)
574 
575  {
576 
577  a[i] = 0.0;
578 
579  b[i] = 0.0;
580 
581  }
582 
583 
584 
585  // Call the routine to calculate the recurrence coefficients
586 
587  RecCoeff(ncoeffs,a,b,alpha,beta);
588 
589 
590 
591  double* a0 = new double[ncoeffs];
592 
593  double* b0 = new double[ncoeffs];
594 
595 
596 
597  chri1(ncoeffs,a,b,a0,b0,end0);
598 
599 
600 
601  double s = b0[0]/fabs(b0[0]);
602 
603  b0[0] = s*b0[0];
604 
605 
606 
607  // Finding the 2*np-1 gauss-kronrod points
608 
609  double* z1 = new double[2*np-1];
610 
611  double* w1 = new double[2*np-1];
612 
613  for(i = 0; i < ncoeffs; i++)
614 
615  {
616 
617  z1[i] = a0[i];
618 
619  w1[i] = b0[i];
620 
621  }
622 
623  JKMatrix(np-1,z1,w1);
624 
625  // Set up the identity matrix
626 
627  double** zmatrix = new double*[2*np-1];
628 
629  for(i = 0; i < 2*np-1; i++)
630 
631  {
632 
633  zmatrix[i] = new double[2*np-1];
634 
635  for(j = 0; j < 2*np-1; j++)
636 
637  {
638 
639  zmatrix[i][j] = 0.0;
640 
641  }
642 
643  }
644 
645  for(i = 0; i < 2*np-1; i++)
646 
647  {
648 
649  zmatrix[i][i] = 1.0;
650 
651  }
652 
653 
654 
655  // Calculate the points and weights
656 
657  TriQL(2*np-1, z1, w1, zmatrix);
658 
659 
660 
661  double sumW1 = 0.0;
662 
663  for(i = 0; i < 2*np-1; i++)
664 
665  {
666 
667  w1[i] = s*w1[i]/(z1[i]-end0);
668 
669  sumW1 += w1[i];
670 
671  }
672 
673 
674 
675  z[0] = end0;
676 
677  w[0] = b[0]- sumW1;
678 
679  for(i = 1; i < kpoints; i++)
680 
681  {
682 
683  z[i] = z1[i-1];
684 
685  w[i] = w1[i-1];
686 
687  }
688 
689 
690 
691 
692 
693  }
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
Definition: Polylib.cpp:2449
Polylib::chri1
void chri1(int, double *, double *, double *, double *, double)
Given a weight function through the first n+1.
Definition: Polylib.cpp:2931
Polylib::JKMatrix
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the.
Definition: Polylib.cpp:2769
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:2565
Generated on Wed Mar 8 2017 12:30:11 for Nektar++ by   doxygen 1.8.8