Nektar++
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 , , 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 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, , and their first derivative, . 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 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 .

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

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  }
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

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  }
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
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

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  }
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
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

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  }
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
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2200
 double Polylib::gammaF ( const double x )

Calculate the Gamma function , , for integer.

values and halves.

Determine the value of using:

where

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.

• Uses the defintion of the Lagrangian interpolant:
%

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  }
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
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:18
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.

• Uses the defintion of the Lagrangian interpolant:
%

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  }
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
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:18
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.

• Uses the defintion of the Lagrangian interpolant:
%

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  }
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
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:18
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.

• Uses the defintion of the Lagrangian interpolant:
%

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  }
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
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:18
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.

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

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

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

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  }
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 at the

np points z.

• To do this we have used the relation

• This formulation is valid for

Definition at line 2146 of file Polylib.cpp.

References jacobfd().

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  }
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, , and their first derivative, .

• This function returns the vectors poly_in and poly_d

containing the value of the order Jacobi polynomial

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]) :

• To calculate the derivative of the polynomial this routine uses

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

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

Definition at line 1946 of file Polylib.cpp.

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.

This routine is only value for

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  }
#define STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
Definition: Polylib.cpp:14
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:18
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

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

where and are the Jacobi (normalized)

orthogonal polynomials ( 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  }
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
Definition: Polylib.cpp:2444
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  }
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  }
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
#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 ,

• 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  }
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
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:40
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 ,

• 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  }
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
Definition: Polylib.cpp:2444
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the.
Definition: Polylib.cpp:2764
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

• 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  }
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:40
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 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 .

• 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  }
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:40
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 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 .

• 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  }
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:40
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 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 ,

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  }
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
Definition: Polylib.cpp:2444
void chri1(int, double *, double *, double *, double *, double)
Given a weight function through the first n+1.
Definition: Polylib.cpp:2926
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the.
Definition: Polylib.cpp:2764
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 ,

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  }
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
Definition: Polylib.cpp:2444
void chri1(int, double *, double *, double *, double *, double)
Given a weight function through the first n+1.
Definition: Polylib.cpp:2926
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the.
Definition: Polylib.cpp:2764
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:2560