Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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...
 

Detailed Description

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

Function Documentation

void Polylib::chri1 ( int  n,
double *  a,
double *  b,
double *  a0,
double *  b0,
double  z 
)

Given a weight function $w(t)$ through the first n+1.

coefficients a and b of its orthogonal polynomials

this routine generates the first n recurrence coefficients for the orthogonal

polynomials relative to the modified weight function $(t-z)w(t)$.

The result will be placed in the array a0 and b0.

Definition at line 2926 of file Polylib.cpp.

Referenced by zwlk(), and zwrk().

2930  {
2931 
2932 
2933 
2934  double q = ceil(3.0*n/2);
2935 
2936  int size = (int)q+1;
2937 
2938  if(size < n+1)
2939 
2940  {
2941 
2942  fprintf(stderr,"input arrays a and b are too short\n");
2943 
2944  }
2945 
2946  double* r = new double[n+1];
2947 
2948  r[0] = z - a[0];
2949 
2950  r[1] = z - a[1] - b[1]/r[0];
2951 
2952  a0[0] = a[1] + r[1] - r[0];
2953 
2954  b0[0] = -r[0]*b[0];
2955 
2956 
2957 
2958  if(n == 1)
2959 
2960  {
2961 
2962  return;
2963 
2964  }
2965 
2966  int k = 0;
2967 
2968  for(k = 1; k < n; k++)
2969 
2970  {
2971 
2972  r[k+1] = z - a[k+1] - b[k+1]/r[k];
2973 
2974  a0[k] = a[k+1] + r[k+1] - r[k];
2975 
2976  b0[k] = b[k] * r[k]/r[k-1];
2977 
2978  }
2979 
2980 
2981 
2982  }
void Polylib::Dgj ( double *  D,
const double *  z,
const int  np,
const double  alpha,
const double  beta 
)

Compute the Derivative Matrix and its transpose associated.

with the Gauss-Jacobi zeros.

  • Compute the derivative matrix, d, and its transpose, dt,

associated with the n_th order Lagrangian interpolants through the

np Gauss-Jacobi points z such that
$ \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) $

Definition at line 932 of file Polylib.cpp.

References jacobd().

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

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

Compute the Derivative Matrix associated with the.

Gauss-Lobatto-Jacobi zeros.

  • Compute the derivative matrix, d, associated with the n_th

order Lagrange interpolants through the np

Gauss-Lobatto-Jacobi points z such that
$ \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) $

Definition at line 1216 of file Polylib.cpp.

References gammaF(), and jacobd().

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

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

Compute the Derivative Matrix and its transpose associated.

with the Gauss-Radau-Jacobi zeros with a zero at z=-1.

  • Compute the derivative matrix, d, associated with the n_th

order Lagrangian interpolants through the np Gauss-Radau-Jacobi

points z such that
$ \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) $

Definition at line 1018 of file Polylib.cpp.

References gammaF(), and jacobd().

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

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

Compute the Derivative Matrix associated with the.

Gauss-Radau-Jacobi zeros with a zero at z=1.

  • Compute the derivative matrix, d, associated with the n_th

order Lagrangian interpolants through the np Gauss-Radau-Jacobi

points z such that
$ \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) $

Definition at line 1116 of file Polylib.cpp.

References gammaF(), and jacobd().

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

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

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

values and halves.

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

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

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

Definition at line 2200 of file Polylib.cpp.

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

2200  {
2201 
2202  double gamma = 1.0;
2203 
2204 
2205 
2206  if (x == -0.5) gamma = -2.0*sqrt(M_PI);
2207 
2208  else if (!x) return gamma;
2209 
2210  else if ((x-(int)x) == 0.5){
2211 
2212  int n = (int) x;
2213 
2214  double tmp = x;
2215 
2216 
2217 
2218  gamma = sqrt(M_PI);
2219 
2220  while(n--){
2221 
2222  tmp -= 1.0;
2223 
2224  gamma *= tmp;
2225 
2226  }
2227 
2228  }
2229 
2230  else if ((x-(int)x) == 0.0){
2231 
2232  int n = (int) x;
2233 
2234  double tmp = x;
2235 
2236 
2237 
2238  while(--n){
2239 
2240  tmp -= 1.0;
2241 
2242  gamma *= tmp;
2243 
2244  }
2245 
2246  }
2247 
2248  else
2249 
2250  fprintf(stderr,"%lf is not of integer or half order\n",x);
2251 
2252  return gamma;
2253 
2254  }
double Polylib::hgj ( const int  i,
const double  z,
const double *  zgj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through.

the np Gauss-Jacobi points zgj at the arbitrary location z.

  • $ -1 \leq z \leq 1 $
  • Uses the defintion of the Lagrangian interpolant:
    %

$ \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{P_{np}^{\alpha,\beta}(z)} {[P_{np}^{\alpha,\beta}(z_j)]^\prime (z-z_j)} & \mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} $

Definition at line 1340 of file Polylib.cpp.

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

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

1344  {
1345 
1346 
1347 
1348  double zi, dz, p, pd, h;
1349 
1350 
1351 
1352  zi = *(zgj+i);
1353 
1354  dz = z - zi;
1355 
1356  if (fabs(dz) < EPS) return 1.0;
1357 
1358 
1359 
1360  jacobd (1, &zi, &pd , np, alpha, beta);
1361 
1362  jacobfd(1, &z , &p, NULL , np, alpha, beta);
1363 
1364  h = p/(pd*dz);
1365 
1366 
1367 
1368  return h;
1369 
1370  }
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:2146
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:18
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:1946
double Polylib::hglj ( const int  i,
const double  z,
const double *  zglj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through the.

np Gauss-Lobatto-Jacobi points zgrj at the arbitrary location

z.

  • $ -1 \leq z \leq 1 $
  • Uses the defintion of the Lagrangian interpolant:
    %

$ \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{(1-z^2) P_{np-2}^{\alpha+1,\beta+1}(z)} {((1-z^2_j) [P_{np-2}^{\alpha+1,\beta+1}(z_j)]^\prime - 2 z_j P_{np-2}^{\alpha+1,\beta+1}(z_j) ) (z-z_j)}&\mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} $

Definition at line 1584 of file Polylib.cpp.

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

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

1588  {
1589 
1590  double one = 1., two = 2.;
1591 
1592  double zi, dz, p, pd, h;
1593 
1594 
1595 
1596  zi = *(zglj+i);
1597 
1598  dz = z - zi;
1599 
1600  if (fabs(dz) < EPS) return 1.0;
1601 
1602 
1603 
1604  jacobfd(1, &zi, &p , NULL, np-2, alpha + one, beta + one);
1605 
1606  // need to use this routine in caes z = -1 or 1
1607 
1608  jacobd (1, &zi, &pd, np-2, alpha + one, beta + one);
1609 
1610  h = (one - zi*zi)*pd - two*zi*p;
1611 
1612  jacobfd(1, &z, &p, NULL, np-2, alpha + one, beta + one);
1613 
1614  h = (one - z*z)*p/(h*dz);
1615 
1616 
1617 
1618  return h;
1619 
1620  }
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:2146
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:18
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:1946
double Polylib::hgrjm ( const int  i,
const double  z,
const double *  zgrj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through the.

np Gauss-Radau-Jacobi points zgrj at the arbitrary location

z. This routine assumes zgrj includes the point -1.

  • $ -1 \leq z \leq 1 $
  • Uses the defintion of the Lagrangian interpolant:
    %

$ \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{(1+z) P_{np-1}^{\alpha,\beta+1}(z)} {((1+z_j) [P_{np-1}^{\alpha,\beta+1}(z_j)]^\prime + P_{np-1}^{\alpha,\beta+1}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} $

Definition at line 1416 of file Polylib.cpp.

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

Referenced by Imgrjm().

1420  {
1421 
1422 
1423 
1424  double zi, dz, p, pd, h;
1425 
1426 
1427 
1428  zi = *(zgrj+i);
1429 
1430  dz = z - zi;
1431 
1432  if (fabs(dz) < EPS) return 1.0;
1433 
1434 
1435 
1436  jacobfd (1, &zi, &p , NULL, np-1, alpha, beta + 1);
1437 
1438  // need to use this routine in caes zi = -1 or 1
1439 
1440  jacobd (1, &zi, &pd, np-1, alpha, beta + 1);
1441 
1442  h = (1.0 + zi)*pd + p;
1443 
1444  jacobfd (1, &z, &p, NULL, np-1, alpha, beta + 1);
1445 
1446  h = (1.0 + z )*p/(h*dz);
1447 
1448 
1449 
1450  return h;
1451 
1452  }
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:2146
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:18
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:1946
double Polylib::hgrjp ( const int  i,
const double  z,
const double *  zgrj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through the.

np Gauss-Radau-Jacobi points zgrj at the arbitrary location

z. This routine assumes zgrj includes the point +1.

  • $ -1 \leq z \leq 1 $
  • Uses the defintion of the Lagrangian interpolant:
    %

$ \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{(1-z) P_{np-1}^{\alpha+1,\beta}(z)} {((1-z_j) [P_{np-1}^{\alpha+1,\beta}(z_j)]^\prime - P_{np-1}^{\alpha+1,\beta}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} $

Definition at line 1500 of file Polylib.cpp.

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

Referenced by Imgrjp().

1504  {
1505 
1506 
1507 
1508  double zi, dz, p, pd, h;
1509 
1510 
1511 
1512  zi = *(zgrj+i);
1513 
1514  dz = z - zi;
1515 
1516  if (fabs(dz) < EPS) return 1.0;
1517 
1518 
1519 
1520  jacobfd (1, &zi, &p , NULL, np-1, alpha+1, beta );
1521 
1522  // need to use this routine in caes z = -1 or 1
1523 
1524  jacobd (1, &zi, &pd, np-1, alpha+1, beta );
1525 
1526  h = (1.0 - zi)*pd - p;
1527 
1528  jacobfd (1, &z, &p, NULL, np-1, alpha+1, beta);
1529 
1530  h = (1.0 - z )*p/(h*dz);
1531 
1532 
1533 
1534  return h;
1535 
1536  }
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:2146
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:18
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:1946
void Polylib::Imgj ( double *  im,
const double *  zgj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Jacobi points to an.

arbitrary distribution at points zm

  • Computes the one-dimensional interpolation matrix, im, to

interpolate a function from at Gauss-Jacobi distribution of nz

zeros zgrj to an arbitrary distribution of mz points zm, i.e.
$ u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) $

Definition at line 1652 of file Polylib.cpp.

References hgj().

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

1654  {
1655 
1656  double zp;
1657 
1658  register int i, j;
1659 
1660 
1661 
1662  for (i = 0; i < nz; ++i) {
1663 
1664  for (j = 0; j < mz; ++j)
1665 
1666  {
1667 
1668  zp = zm[j];
1669 
1670  im [i*mz+j] = hgj(i, zp, zgj, nz, alpha, beta);
1671 
1672  }
1673 
1674  }
1675 
1676 
1677 
1678  return;
1679 
1680  }
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:1340
void Polylib::Imglj ( double *  im,
const double *  zglj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Lobatto-Jacobi points.

to an arbitrary distrubtion at points zm

  • Computes the one-dimensional interpolation matrix, im, to

interpolate a function from at Gauss-Lobatto-Jacobi distribution of

nz zeros zgrj (where zgrj[0]=-1) to an arbitrary

distribution of mz points zm, i.e.


$ u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) $

Definition at line 1832 of file Polylib.cpp.

References hglj().

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

1836  {
1837 
1838  double zp;
1839 
1840  register int i, j;
1841 
1842 
1843 
1844  for (i = 0; i < nz; i++) {
1845 
1846  for (j = 0; j < mz; j++)
1847 
1848  {
1849 
1850  zp = zm[j];
1851 
1852  im[i*mz+j] = hglj(i, zp, zglj, nz, alpha, beta);
1853 
1854  }
1855 
1856  }
1857 
1858 
1859 
1860  return;
1861 
1862  }
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:1584
void Polylib::Imgrjm ( double *  im,
const double *  zgrj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Radau-Jacobi points.

(including z=-1) to an arbitrary distrubtion at points zm

  • Computes the one-dimensional interpolation matrix, im, to

interpolate a function from at Gauss-Radau-Jacobi distribution of

nz zeros zgrj (where zgrj[0]=-1) to an arbitrary

distribution of mz points zm, i.e.


$ u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) $

Definition at line 1710 of file Polylib.cpp.

References hgrjm().

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

1714  {
1715 
1716  double zp;
1717 
1718  register int i, j;
1719 
1720 
1721 
1722  for (i = 0; i < nz; i++) {
1723 
1724  for (j = 0; j < mz; j++)
1725 
1726  {
1727 
1728  zp = zm[j];
1729 
1730  im [i*mz+j] = hgrjm(i, zp, zgrj, nz, alpha, beta);
1731 
1732  }
1733 
1734  }
1735 
1736 
1737 
1738  return;
1739 
1740  }
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:1416
void Polylib::Imgrjp ( double *  im,
const double *  zgrj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Radau-Jacobi points.

(including z=1) to an arbitrary distrubtion at points zm

  • Computes the one-dimensional interpolation matrix, im, to

interpolate a function from at Gauss-Radau-Jacobi distribution of

nz zeros zgrj (where zgrj[nz-1]=1) to an arbitrary

distribution of mz points zm, i.e.


$ u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) $

Definition at line 1770 of file Polylib.cpp.

References hgrjp().

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

1774  {
1775 
1776  double zp;
1777 
1778  register int i, j;
1779 
1780 
1781 
1782  for (i = 0; i < nz; i++) {
1783 
1784  for (j = 0; j < mz; j++)
1785 
1786  {
1787 
1788  zp = zm[j];
1789 
1790  im [i*mz+j] = hgrjp(i, zp, zgrj, nz, alpha, beta);
1791 
1792  }
1793 
1794  }
1795 
1796 
1797 
1798  return;
1799 
1800  }
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:1500
void Polylib::jacobd ( const int  np,
const double *  z,
double *  polyd,
const int  n,
const double  alpha,
const double  beta 
)

Calculate the derivative of Jacobi polynomials.

  • Generates a vector poly of values of the derivative of the

n th order Jacobi polynomial $ P^(\alpha,\beta)_n(z)$ at the

np points z.

  • To do this we have used the relation


$ \frac{d}{dz} P^{\alpha,\beta}_n(z) = \frac{1}{2} (\alpha + \beta + n + 1) P^{\alpha,\beta}_n(z) $

  • This formulation is valid for $ -1 \leq z \leq 1 $

Definition at line 2146 of file Polylib.cpp.

References jacobfd().

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

2150  {
2151 
2152  register int i;
2153 
2154  double one = 1.0;
2155 
2156  if(n == 0)
2157 
2158  for(i = 0; i < np; ++i) polyd[i] = 0.0;
2159 
2160  else{
2161 
2162  //jacobf(np,z,polyd,n-1,alpha+one,beta+one);
2163 
2164  jacobfd(np,z,polyd,NULL,n-1,alpha+one,beta+one);
2165 
2166  for(i = 0; i < np; ++i) polyd[i] *= 0.5*(alpha + beta + (double)n + one);
2167 
2168  }
2169 
2170  return;
2171 
2172  }
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:1946
void Polylib::jacobfd ( const int  np,
const double *  z,
double *  poly_in,
double *  polyd,
const int  n,
const double  alpha,
const double  beta 
)

Routine to calculate Jacobi polynomials, $ P^{\alpha,\beta}_n(z) $, and their first derivative, $ \frac{d}{dz} P^{\alpha,\beta}_n(z) $.

  • This function returns the vectors poly_in and poly_d

containing the value of the $ n^th $ order Jacobi polynomial

$ P^{\alpha,\beta}_n(z) \alpha > -1, \beta > -1 $ and its

derivative at the np points in z[i]

  • If poly_in = NULL then only calculate derivatice
  • If polyd = NULL then only calculate polynomial
  • To calculate the polynomial this routine uses the recursion

relationship (see appendix A ref [4]) :

$ \begin{array}{rcl} P^{\alpha,\beta}_0(z) &=& 1 \\ P^{\alpha,\beta}_1(z) &=& \frac{1}{2} [ \alpha-\beta+(\alpha+\beta+2)z] \\ a^1_n P^{\alpha,\beta}_{n+1}(z) &=& (a^2_n + a^3_n z) P^{\alpha,\beta}_n(z) - a^4_n P^{\alpha,\beta}_{n-1}(z) \\ a^1_n &=& 2(n+1)(n+\alpha + \beta + 1)(2n + \alpha + \beta) \\ a^2_n &=& (2n + \alpha + \beta + 1)(\alpha^2 - \beta^2) \\ a^3_n &=& (2n + \alpha + \beta)(2n + \alpha + \beta + 1) (2n + \alpha + \beta + 2) \\ a^4_n &=& 2(n+\alpha)(n+\beta)(2n + \alpha + \beta + 2) \end{array} $

  • To calculate the derivative of the polynomial this routine uses

the relationship (see appendix A ref [4]) :

$ \begin{array}{rcl} b^1_n(z)\frac{d}{dz} P^{\alpha,\beta}_n(z)&=&b^2_n(z)P^{\alpha,\beta}_n(z) + b^3_n(z) P^{\alpha,\beta}_{n-1}(z) \hspace{2.2cm} \\ b^1_n(z) &=& (2n+\alpha + \beta)(1-z^2) \\ b^2_n(z) &=& n[\alpha - \beta - (2n+\alpha + \beta)z]\\ b^3_n(z) &=& 2(n+\alpha)(n+\beta) \end{array} $

  • Note the derivative from this routine is only valid for -1 < z < 1.

Definition at line 1946 of file Polylib.cpp.

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

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

Calculate the n zeros, z, of the Jacobi polynomial, i.e.

$ P_n^{\alpha,\beta}(z) = 0 $

This routine is only value for $( \alpha > -1, \beta > -1)$

and uses polynomial deflation in a Newton iteration

Definition at line 2274 of file Polylib.cpp.

References EPS, jacobfd(), and STOP.

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

Zero and Weight determination through the eigenvalues and eigenvectors of a tridiagonal.

matrix from the three term recurrence relationship.

Set up a symmetric tridiagonal matrix

$ \left [ \begin{array}{ccccc} a[0] & b[0] & & & \\ b[0] & a[1] & b[1] & & \\ 0 & \ddots & \ddots & \ddots & \\ & & \ddots & \ddots & b[n-2] \\ & & & b[n-2] & a[n-1] \end{array} \right ] $

Where the coefficients a[n], b[n] come from the recurrence relation

$ b_j p_j(z) = (z - a_j ) p_{j-1}(z) - b_{j-1} p_{j-2}(z) $

where $ j=n+1$ and $p_j(z)$ are the Jacobi (normalized)

orthogonal polynomials $ \alpha,\beta > -1$( integer values and

halves). Since the polynomials are orthonormalized, the tridiagonal

matrix is guaranteed to be symmetric. The eigenvalues of this

matrix are the zeros of the Jacobi polynomial.

Definition at line 2384 of file Polylib.cpp.

References RecCoeff(), and TriQL().

2386  {
2387 
2388 
2389 
2390  int i,j;
2391 
2392  RecCoeff(n,a,b,alpha,beta);
2393 
2394 
2395 
2396  double **z = new double*[n];
2397 
2398  for(i = 0; i < n; i++)
2399 
2400  {
2401 
2402  z[i] = new double[n];
2403 
2404  for(j = 0; j < n; j++)
2405 
2406  {
2407 
2408  z[i][j] = 0.0;
2409 
2410  }
2411 
2412  }
2413 
2414  for(i = 0; i < n; i++)
2415 
2416  {
2417 
2418  z[i][i] = 1.0;
2419 
2420  }
2421 
2422 
2423 
2424  // find eigenvalues and eigenvectors
2425 
2426  TriQL(n, a, b,z);
2427 
2428 
2429 
2430  return;
2431 
2432  }
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
Definition: Polylib.cpp:2444
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:2560
void Polylib::JKMatrix ( int  n,
double *  a,
double *  b 
)

Calcualtes the Jacobi-kronrod matrix by determining the.

a and coefficients.

The first 3n+1 coefficients are already known

For more information refer to:

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

Definition at line 2764 of file Polylib.cpp.

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

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

The routine finds the recurrence coefficients a and.

b of the orthogonal polynomials

Definition at line 2444 of file Polylib.cpp.

References gammaF().

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

2446  {
2447 
2448 
2449 
2450  int i;
2451 
2452  double apb, apbi,a2b2;
2453 
2454 
2455 
2456  if(!n)
2457 
2458  return;
2459 
2460 
2461 
2462  // generate normalised terms
2463 
2464  apb = alpha + beta;
2465 
2466  apbi = 2.0 + apb;
2467 
2468 
2469 
2470  b[0] = pow(2.0,apb+1.0)*gammaF(alpha+1.0)*gammaF(beta+1.0)/gammaF(apbi); //MuZero
2471 
2472  a[0] = (beta-alpha)/apbi;
2473 
2474  b[1] = (4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));
2475 
2476 
2477 
2478  a2b2 = beta*beta-alpha*alpha;
2479 
2480 
2481 
2482  for(i = 1; i < n-1; i++){
2483 
2484  apbi = 2.0*(i+1) + apb;
2485 
2486  a[i] = a2b2/((apbi-2.0)*apbi);
2487 
2488  b[i+1] = (4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
2489 
2490  ((apbi*apbi-1)*apbi*apbi));
2491 
2492  }
2493 
2494 
2495 
2496  apbi = 2.0*n + apb;
2497 
2498  a[n-1] = a2b2/((apbi-2.0)*apbi);
2499 
2500 
2501 
2502  }
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2200
static void Polylib::TriQL ( const int  n,
double *  d,
double *  e,
double **  z 
)
static

QL algorithm for symmetric tridiagonal matrix.

This subroutine is a translation of an algol procedure,

num. math. 12, 377-383(1968) by martin and wilkinson, as modified

in num. math. 15, 450(1970) by dubrulle. Handbook for

auto. comp., vol.ii-linear algebra, 241-248(1971). This is a

modified version from numerical recipes.

This subroutine finds the eigenvalues and first components of the

eigenvectors of a symmetric tridiagonal matrix by the implicit QL

method.

on input:

  • n is the order of the matrix;
  • d contains the diagonal elements of the input matrix;
  • e contains the subdiagonal elements of the input matrix

in its first n-2 positions.

  • z is the n by n identity matrix

on output:

  • d contains the eigenvalues in ascending order.

  • e contains the weight values - modifications of the first component

    of normalised eigenvectors

Definition at line 2560 of file Polylib.cpp.

References sign, and STOP.

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

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

Gauss-Jacobi zeros and weights.

  • Generate np Gauss Jacobi zeros, z, and weights,w,

associated with the Jacobi polynomial $ P^{\alpha,\beta}_{np}(z) $,

  • Exact for polynomials of order 2np-1 or less

Definition at line 102 of file Polylib.cpp.

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

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

106  {
107 
108  register int i;
109 
110  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
111 
112 
113 
114  jacobz (np,z,alpha,beta);
115 
116  jacobd (np,z,w,np,alpha,beta);
117 
118 
119 
120  fac = pow(two,apb + one)*gammaF(alpha + np + one)*gammaF(beta + np + one);
121 
122  fac /= gammaF(np + one)*gammaF(apb + np + one);
123 
124 
125 
126  for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]*(one-z[i]*z[i]));
127 
128 
129 
130  return;
131 
132  }
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:2146
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:40
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2200
void Polylib::zwgk ( double *  z,
double *  w,
const int  npt,
const double  alpha,
const double  beta 
)

Gauss-Kronrod-Jacobi zeros and weights.

  • Generate npt=2*np+1 Gauss-Kronrod Jacobi zeros, z, and weights,w,

associated with the Jacobi polynomial $ P^{\alpha,\beta}_{np}(z) $,

  • Exact for polynomials of order 3np+1 or less

Definition at line 380 of file Polylib.cpp.

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

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

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

Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1.

  • Generate np Gauss-Lobatto-Jacobi points, z, and weights, w,

associated with polynomial $ (1-z)(1+z) P^{\alpha+1,\beta+1}_{np-2}(z) $

  • Exact for polynomials of order 2np-3 or less

Definition at line 306 of file Polylib.cpp.

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

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

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

Gauss-Radau-Jacobi zeros and weights with end point at z=-1.

  • Generate np Gauss-Radau-Jacobi zeros, z, and weights,w,

associated with the polynomial $(1+z) P^{\alpha,\beta+1}_{np-1}(z) $.

  • Exact for polynomials of order 2np-2 or less

Definition at line 158 of file Polylib.cpp.

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

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

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

Gauss-Radau-Jacobi zeros and weights with end point at z=1.

  • Generate np Gauss-Radau-Jacobi zeros, z, and weights,w,

associated with the polynomial $(1-z) P^{\alpha+1,\beta}_{np-1}(z) $.

  • Exact for polynomials of order 2np-2 or less

Definition at line 234 of file Polylib.cpp.

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

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

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

Gauss-Lobatto-Kronrod-Jacobi zeros and weights.

  • Generate npt=2*np-1 Lobatto-Kronrod Jacobi zeros, z, and weights,w,

associated with the Jacobi polynomial $ P^{\alpha,\beta}_{np}(z) $,

Definition at line 706 of file Polylib.cpp.

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

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

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

Gauss-Radau-Kronrod-Jacobi zeros and weights.

  • Generate npt=2*np Radau-Kronrod Jacobi zeros, z, and weights,w,

associated with the Jacobi polynomial $ P^{\alpha,\beta}_{np}(z) $,

Definition at line 514 of file Polylib.cpp.

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

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

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