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

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

Functions

double optdiff (double xl, double xr)
 The following function is used to circumvent/reduce "Subtractive Cancellation" The expression 1/dz is replaced by optinvsub(.,.) Added on 26 April 2017. More...
 
double laginterp (double z, int j, const double *zj, int np)
 
static void Jacobz (const int n, double *z, const double alpha, const double beta)
 Calculate the n zeros, z, of the Jacobi polynomial, i.e. More...
 
static void TriQL (const int n, double *d, double *e, double **z)
 QL algorithm for symmetric tridiagonal matrix. More...
 
double gammaF (const double x)
 Calculate the Gamma function , $ \Gamma(n)$, for integer. More...
 
static void RecCoeff (const int n, double *a, double *b, const double alpha, const double beta)
 The routine finds the recurrence coefficients a and. More...
 
void JKMatrix (int n, double *a, double *b)
 Calcualtes the Jacobi-kronrod matrix by determining the. More...
 
void chri1 (int n, double *a, double *b, double *a0, double *b0, double z)
 Given a weight function $w(t)$ through the first n+1. More...
 
void zwgj (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Jacobi zeros and weights. More...
 
void zwgrjm (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Radau-Jacobi zeros and weights with end point at z=-1. More...
 
void zwgrjp (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Radau-Jacobi zeros and weights with end point at z=1. More...
 
void zwglj (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1. More...
 
void zwgk (double *z, double *w, const int npt, const double alpha, const double beta)
 Gauss-Kronrod-Jacobi zeros and weights. More...
 
void zwrk (double *z, double *w, const int npt, const double alpha, const double beta)
 Gauss-Radau-Kronrod-Jacobi zeros and weights. More...
 
void zwlk (double *z, double *w, const int npt, const double alpha, const double beta)
 Gauss-Lobatto-Kronrod-Jacobi zeros and weights. More...
 
void Dgj (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix and its transpose associated. More...
 
void Dgrjm (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix and its transpose associated. More...
 
void Dgrjp (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix associated with the. More...
 
void Dglj (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix associated with the. More...
 
double hgj (const int i, const double z, const double *zgj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through. More...
 
double hgrjm (const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through the. More...
 
double hgrjp (const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through the. More...
 
double hglj (const int i, const double z, const double *zglj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through the. More...
 
void Imgj (double *im, const double *zgj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Jacobi points to an. More...
 
void Imgrjm (double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Radau-Jacobi points. More...
 
void Imgrjp (double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Radau-Jacobi points. More...
 
void Imglj (double *im, const double *zglj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Lobatto-Jacobi points. More...
 
void jacobfd (const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
 Routine to calculate Jacobi polynomials, $ P^{\alpha,\beta}_n(z) $, and their first derivative, $ \frac{d}{dz} P^{\alpha,\beta}_n(z) $. More...
 
void jacobd (const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
 Calculate the derivative of Jacobi polynomials. More...
 
void JacZeros (const int n, double *a, double *b, const double alpha, const double beta)
 Zero and Weight determination through the eigenvalues and eigenvectors of a tridiagonal. More...
 
std::complex< Nektar::NekDouble > ImagBesselComp (int n, std::complex< Nektar::NekDouble > y)
 Calcualte the bessel function of the first kind with complex double input y. Taken from Numerical Recipies in C. More...
 

Detailed Description

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

Function Documentation

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

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

coefficients a and b of its orthogonal polynomials

this routine generates the first n recurrence coefficients for the orthogonal

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

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

Definition at line 2900 of file Polylib.cpp.

Referenced by zwlk(), and zwrk().

2904  {
2905 
2906 
2907 
2908  double q = ceil(3.0*n/2);
2909 
2910  int size = (int)q+1;
2911 
2912  if(size < n+1)
2913 
2914  {
2915 
2916  fprintf(stderr,"input arrays a and b are too short\n");
2917 
2918  }
2919 
2920  double* r = new double[n+1];
2921 
2922  r[0] = z - a[0];
2923 
2924  r[1] = z - a[1] - b[1]/r[0];
2925 
2926  a0[0] = a[1] + r[1] - r[0];
2927 
2928  b0[0] = -r[0]*b[0];
2929 
2930 
2931 
2932  if(n == 1)
2933 
2934  {
2935 
2936  return;
2937 
2938  }
2939 
2940  int k = 0;
2941 
2942  for(k = 1; k < n; k++)
2943 
2944  {
2945 
2946  r[k+1] = z - a[k+1] - b[k+1]/r[k];
2947 
2948  a0[k] = a[k+1] + r[k+1] - r[k];
2949 
2950  b0[k] = b[k] * r[k]/r[k-1];
2951 
2952  }
2953 
2954 
2955  }
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 972 of file Polylib.cpp.

References jacobd().

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

976  {
977 
978 
979 
980  double one = 1.0, two = 2.0;
981 
982 
983 
984  if (np <= 0){
985 
986  D[0] = 0.0;
987 
988  }
989 
990  else{
991 
992  register int i,j;
993 
994  double *pd;
995 
996 
997 
998  pd = (double *)malloc(np*sizeof(double));
999 
1000  jacobd(np,z,pd,np,alpha,beta);
1001 
1002 
1003 
1004  for (i = 0; i < np; i++){
1005 
1006  for (j = 0; j < np; j++){
1007 
1008 
1009 
1010  if (i != j)
1011 
1012  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1013 
1014  else
1015 
1016  D[i*np+j] = (alpha - beta + (alpha + beta + two)*z[j])/
1017 
1018  (two*(one - z[j]*z[j]));
1019 
1020  }
1021 
1022  }
1023 
1024  free(pd);
1025 
1026  }
1027 
1028  return;
1029 
1030  }
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:2120
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 1260 of file Polylib.cpp.

References gammaF(), and jacobd().

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

1264  {
1265 
1266 
1267 
1268  if (np <= 0){
1269 
1270  D[0] = 0.0;
1271 
1272  }
1273 
1274  else{
1275 
1276  register int i, j;
1277 
1278  double one = 1.0, two = 2.0;
1279 
1280  double *pd;
1281 
1282 
1283 
1284  pd = (double *)malloc(np*sizeof(double));
1285 
1286 
1287 
1288  pd[0] = two*pow(-one,np)*gammaF(np + beta);
1289 
1290  pd[0] /= gammaF(np - one)*gammaF(beta + two);
1291 
1292  jacobd(np-2,z+1,pd+1,np-2,alpha+1,beta+1);
1293 
1294  for(i = 1; i < np-1; ++i) pd[i] *= (one-z[i]*z[i]);
1295 
1296  pd[np-1] = -two*gammaF(np + alpha);
1297 
1298  pd[np-1] /= gammaF(np - one)*gammaF(alpha + two);
1299 
1300 
1301 
1302  for (i = 0; i < np; i++) {
1303 
1304  for (j = 0; j < np; j++){
1305 
1306  if (i != j)
1307 
1308  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1309 
1310  else {
1311 
1312  if (j == 0)
1313 
1314  D[i*np+j] = (alpha - (np-1)*(np + alpha + beta))/(two*(beta+ two));
1315 
1316  else if (j == np-1)
1317 
1318  D[i*np+j] =-(beta - (np-1)*(np + alpha + beta))/(two*(alpha+ two));
1319 
1320  else
1321 
1322  D[i*np+j] = (alpha - beta + (alpha + beta)*z[j])/
1323 
1324  (two*(one - z[j]*z[j]));
1325 
1326  }
1327 
1328  }
1329 
1330  }
1331 
1332  free(pd);
1333 
1334  }
1335 
1336 
1337 
1338  return;
1339 
1340  }
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:2120
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2174
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 1058 of file Polylib.cpp.

References gammaF(), and jacobd().

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

1062  {
1063 
1064 
1065 
1066  if (np <= 0){
1067 
1068  D[0] = 0.0;
1069 
1070  }
1071 
1072  else{
1073 
1074  register int i, j;
1075 
1076  double one = 1.0, two = 2.0;
1077 
1078  double *pd;
1079 
1080 
1081 
1082  pd = (double *)malloc(np*sizeof(double));
1083 
1084 
1085 
1086  pd[0] = pow(-one,np-1)*gammaF(np+beta+one);
1087 
1088  pd[0] /= gammaF(np)*gammaF(beta+two);
1089 
1090  jacobd(np-1,z+1,pd+1,np-1,alpha,beta+1);
1091 
1092  for(i = 1; i < np; ++i) pd[i] *= (1+z[i]);
1093 
1094 
1095 
1096  for (i = 0; i < np; i++) {
1097 
1098  for (j = 0; j < np; j++){
1099 
1100  if (i != j)
1101 
1102  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1103 
1104  else {
1105 
1106  if(j == 0)
1107 
1108  D[i*np+j] = -(np + alpha + beta + one)*(np - one)/
1109 
1110  (two*(beta + two));
1111 
1112  else
1113 
1114  D[i*np+j] = (alpha - beta + one + (alpha + beta + one)*z[j])/
1115 
1116  (two*(one - z[j]*z[j]));
1117 
1118  }
1119 
1120  }
1121 
1122  }
1123 
1124  free(pd);
1125 
1126  }
1127 
1128 
1129 
1130  return;
1131 
1132  }
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:2120
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2174
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 1158 of file Polylib.cpp.

References gammaF(), and jacobd().

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

1162  {
1163 
1164 
1165 
1166  if (np <= 0){
1167 
1168  D[0] = 0.0;
1169 
1170  }
1171 
1172  else{
1173 
1174  register int i, j;
1175 
1176  double one = 1.0, two = 2.0;
1177 
1178  double *pd;
1179 
1180 
1181 
1182  pd = (double *)malloc(np*sizeof(double));
1183 
1184 
1185 
1186 
1187 
1188  jacobd(np-1,z,pd,np-1,alpha+1,beta);
1189 
1190  for(i = 0; i < np-1; ++i) pd[i] *= (1-z[i]);
1191 
1192  pd[np-1] = -gammaF(np+alpha+one);
1193 
1194  pd[np-1] /= gammaF(np)*gammaF(alpha+two);
1195 
1196 
1197 
1198  for (i = 0; i < np; i++) {
1199 
1200  for (j = 0; j < np; j++){
1201 
1202  if (i != j)
1203 
1204  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1205 
1206  else {
1207 
1208  if(j == np-1)
1209 
1210  D[i*np+j] = (np + alpha + beta + one)*(np - one)/
1211 
1212  (two*(alpha + two));
1213 
1214  else
1215 
1216  D[i*np+j] = (alpha - beta - one + (alpha + beta + one)*z[j])/
1217 
1218  (two*(one - z[j]*z[j]));
1219 
1220  }
1221 
1222  }
1223 
1224  }
1225 
1226  free(pd);
1227 
1228  }
1229 
1230 
1231 
1232  return;
1233 
1234  }
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:2120
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2174
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 2174 of file Polylib.cpp.

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

2174  {
2175 
2176  double gamma = 1.0;
2177 
2178 
2179 
2180  if (x == -0.5) gamma = -2.0*sqrt(M_PI);
2181 
2182  else if (!x) return gamma;
2183 
2184  else if ((x-(int)x) == 0.5){
2185 
2186  int n = (int) x;
2187 
2188  double tmp = x;
2189 
2190 
2191 
2192  gamma = sqrt(M_PI);
2193 
2194  while(n--){
2195 
2196  tmp -= 1.0;
2197 
2198  gamma *= tmp;
2199 
2200  }
2201 
2202  }
2203 
2204  else if ((x-(int)x) == 0.0){
2205 
2206  int n = (int) x;
2207 
2208  double tmp = x;
2209 
2210 
2211 
2212  while(--n){
2213 
2214  tmp -= 1.0;
2215 
2216  gamma *= tmp;
2217 
2218  }
2219 
2220  }
2221 
2222  else
2223 
2224  fprintf(stderr,"%lf is not of integer or half order\n",x);
2225 
2226  return gamma;
2227 
2228  }
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 1386 of file Polylib.cpp.

References EPS, and laginterp().

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

1390  {
1391 
1392  double zi, dz;
1393 
1394  zi = *(zgj+i);
1395 
1396  dz = z-zi;
1397 
1398  if (fabs(dz) < EPS) return 1.0;
1399 
1400  return laginterp(z, i, zgj, np);
1401 
1402  }
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:23
Polylib::laginterp
double laginterp(double z, int j, const double *zj, int np)
Definition: Polylib.cpp:57
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 1578 of file Polylib.cpp.

References EPS, and laginterp().

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

1582  {
1583 
1584  double zi, dz;
1585 
1586  zi = *(zglj+i);
1587 
1588  dz = z-zi;
1589 
1590  if (fabs(dz) < EPS) return 1.0;
1591 
1592  return laginterp(z, i, zglj, np);
1593 
1594  }
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:23
Polylib::laginterp
double laginterp(double z, int j, const double *zj, int np)
Definition: Polylib.cpp:57
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 1448 of file Polylib.cpp.

References EPS, and laginterp().

Referenced by Imgrjm().

1452  {
1453 
1454 
1455  double zi, dz;
1456 
1457  zi = *(zgrj+i);
1458 
1459  dz = z-zi;
1460 
1461  if (fabs(dz) < EPS) return 1.0;
1462 
1463  return laginterp(z, i, zgrj, np);
1464 
1465  }
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:23
Polylib::laginterp
double laginterp(double z, int j, const double *zj, int np)
Definition: Polylib.cpp:57
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 1513 of file Polylib.cpp.

References EPS, and laginterp().

Referenced by Imgrjp().

1517  {
1518 
1519 
1520  double zi, dz;
1521 
1522  zi = *(zgrj+i);
1523 
1524  dz = z-zi;
1525 
1526  if (fabs(dz) < EPS) return 1.0;
1527 
1528  return laginterp(z, i, zgrj, np);
1529 
1530  }
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:23
Polylib::laginterp
double laginterp(double z, int j, const double *zj, int np)
Definition: Polylib.cpp:57
std::complex< Nektar::NekDouble > Polylib::ImagBesselComp ( int  n,
std::complex< Nektar::NekDouble >  y 
)

Calcualte the bessel function of the first kind with complex double input y. Taken from Numerical Recipies in C.

Returns a complex double

Definition at line 2968 of file Polylib.cpp.

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

2969  {
2970  std::complex<Nektar::NekDouble> z (1.0,0.0);
2971  std::complex<Nektar::NekDouble> zbes (1.0,0.0);
2972  std::complex<Nektar::NekDouble> zarg;
2973  Nektar::NekDouble tol = 1e-15;
2974  int maxit = 10000;
2975  int i = 1;
2976 
2977  zarg = -0.25*y*y;
2978 
2979  while (abs(z) > tol && i <= maxit){
2980  z = z*(1.0/i/(i+n)*zarg);
2981  if (abs(z) <= tol) break;
2982  zbes = zbes + z;
2983  i++;
2984  }
2985  zarg = 0.5*y;
2986  for (i=1;i<=n;i++){
2987  zbes = zbes*zarg;
2988  }
2989  return zbes;
2990 
2991  }
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
void Polylib::Imgj ( double *  im,
const double *  zgj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Jacobi points to an.

arbitrary distribution at points zm

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

interpolate a function from at Gauss-Jacobi distribution of nz

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

Definition at line 1626 of file Polylib.cpp.

References hgj().

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

1628  {
1629 
1630  double zp;
1631 
1632  register int i, j;
1633 
1634 
1635 
1636  for (i = 0; i < nz; ++i) {
1637 
1638  for (j = 0; j < mz; ++j)
1639 
1640  {
1641 
1642  zp = zm[j];
1643 
1644  im [i*mz+j] = hgj(i, zp, zgj, nz, alpha, beta);
1645 
1646  }
1647 
1648  }
1649 
1650 
1651 
1652  return;
1653 
1654  }
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:1386
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 1806 of file Polylib.cpp.

References hglj().

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

1810  {
1811 
1812  double zp;
1813 
1814  register int i, j;
1815 
1816 
1817 
1818  for (i = 0; i < nz; i++) {
1819 
1820  for (j = 0; j < mz; j++)
1821 
1822  {
1823 
1824  zp = zm[j];
1825 
1826  im[i*mz+j] = hglj(i, zp, zglj, nz, alpha, beta);
1827 
1828  }
1829 
1830  }
1831 
1832 
1833 
1834  return;
1835 
1836  }
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:1578
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 1684 of file Polylib.cpp.

References hgrjm().

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

1688  {
1689 
1690  double zp;
1691 
1692  register int i, j;
1693 
1694 
1695 
1696  for (i = 0; i < nz; i++) {
1697 
1698  for (j = 0; j < mz; j++)
1699 
1700  {
1701 
1702  zp = zm[j];
1703 
1704  im [i*mz+j] = hgrjm(i, zp, zgrj, nz, alpha, beta);
1705 
1706  }
1707 
1708  }
1709 
1710 
1711 
1712  return;
1713 
1714  }
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:1448
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 1744 of file Polylib.cpp.

References hgrjp().

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

1748  {
1749 
1750  double zp;
1751 
1752  register int i, j;
1753 
1754 
1755 
1756  for (i = 0; i < nz; i++) {
1757 
1758  for (j = 0; j < mz; j++)
1759 
1760  {
1761 
1762  zp = zm[j];
1763 
1764  im [i*mz+j] = hgrjp(i, zp, zgrj, nz, alpha, beta);
1765 
1766  }
1767 
1768  }
1769 
1770 
1771 
1772  return;
1773 
1774  }
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:1513
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 2120 of file Polylib.cpp.

References jacobfd().

Referenced by Dgj(), Dglj(), Dgrjm(), Dgrjp(), Nektar::LibUtilities::NodalUtilTriangle::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilTetrahedron::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilPrism::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilQuad::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilHex::v_OrthoBasisDeriv(), Nektar::SolverUtils::AdvectionFR::v_SetupCFunctions(), Nektar::SolverUtils::DiffusionLFR::v_SetupCFunctions(), Nektar::SolverUtils::DiffusionLFRNS::v_SetupCFunctions(), and zwgj().

2124  {
2125 
2126  register int i;
2127 
2128  double one = 1.0;
2129 
2130  if(n == 0)
2131 
2132  for(i = 0; i < np; ++i) polyd[i] = 0.0;
2133 
2134  else{
2135 
2136  //jacobf(np,z,polyd,n-1,alpha+one,beta+one);
2137 
2138  jacobfd(np,z,polyd,NULL,n-1,alpha+one,beta+one);
2139 
2140  for(i = 0; i < np; ++i) polyd[i] *= 0.5*(alpha + beta + (double)n + one);
2141 
2142  }
2143 
2144  return;
2145 
2146  }
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:1920
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 1920 of file Polylib.cpp.

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

1922  {
1923 
1924  register int i;
1925 
1926  double zero = 0.0, one = 1.0, two = 2.0;
1927 
1928 
1929 
1930  if(!np)
1931 
1932  return;
1933 
1934 
1935 
1936  if(n == 0){
1937 
1938  if(poly_in)
1939 
1940  for(i = 0; i < np; ++i)
1941 
1942  poly_in[i] = one;
1943 
1944  if(polyd)
1945 
1946  for(i = 0; i < np; ++i)
1947 
1948  polyd[i] = zero;
1949 
1950  }
1951 
1952  else if (n == 1){
1953 
1954  if(poly_in)
1955 
1956  for(i = 0; i < np; ++i)
1957 
1958  poly_in[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
1959 
1960  if(polyd)
1961 
1962  for(i = 0; i < np; ++i)
1963 
1964  polyd[i] = 0.5*(alpha + beta + two);
1965 
1966  }
1967 
1968  else{
1969 
1970  register int k;
1971 
1972  double a1,a2,a3,a4;
1973 
1974  double two = 2.0, apb = alpha + beta;
1975 
1976  double *poly, *polyn1,*polyn2;
1977 
1978 
1979 
1980  if(poly_in){ // switch for case of no poynomial function return
1981 
1982  polyn1 = (double *)malloc(2*np*sizeof(double));
1983 
1984  polyn2 = polyn1+np;
1985 
1986  poly = poly_in;
1987 
1988  }
1989 
1990  else{
1991 
1992  polyn1 = (double *)malloc(3*np*sizeof(double));
1993 
1994  polyn2 = polyn1+np;
1995 
1996  poly = polyn2+np;
1997 
1998  }
1999 
2000 
2001 
2002  for(i = 0; i < np; ++i){
2003 
2004  polyn2[i] = one;
2005 
2006  polyn1[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
2007 
2008  }
2009 
2010 
2011 
2012  for(k = 2; k <= n; ++k){
2013 
2014  a1 = two*k*(k + apb)*(two*k + apb - two);
2015 
2016  a2 = (two*k + apb - one)*(alpha*alpha - beta*beta);
2017 
2018  a3 = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb);
2019 
2020  a4 = two*(k + alpha - one)*(k + beta - one)*(two*k + apb);
2021 
2022 
2023 
2024  a2 /= a1;
2025 
2026  a3 /= a1;
2027 
2028  a4 /= a1;
2029 
2030 
2031 
2032  for(i = 0; i < np; ++i){
2033 
2034  poly [i] = (a2 + a3*z[i])*polyn1[i] - a4*polyn2[i];
2035 
2036  polyn2[i] = polyn1[i];
2037 
2038  polyn1[i] = poly [i];
2039 
2040  }
2041 
2042  }
2043 
2044 
2045 
2046  if(polyd){
2047 
2048  a1 = n*(alpha - beta);
2049 
2050  a2 = n*(two*n + alpha + beta);
2051 
2052  a3 = two*(n + alpha)*(n + beta);
2053 
2054  a4 = (two*n + alpha + beta);
2055 
2056  a1 /= a4; a2 /= a4; a3 /= a4;
2057 
2058 
2059 
2060  // note polyn2 points to polyn1 at end of poly iterations
2061 
2062  for(i = 0; i < np; ++i){
2063 
2064  polyd[i] = (a1- a2*z[i])*poly[i] + a3*polyn2[i];
2065 
2066  polyd[i] /= (one - z[i]*z[i]);
2067 
2068  }
2069 
2070  }
2071 
2072 
2073 
2074  free(polyn1);
2075 
2076  }
2077 
2078 
2079 
2080  return;
2081 
2082  }
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 2248 of file Polylib.cpp.

References EPS, jacobfd(), and STOP.

2250  {
2251 
2252  register int i,j,k;
2253 
2254  double dth = M_PI/(2.0*(double)n);
2255 
2256  double poly,pder,rlast=0.0;
2257 
2258  double sum,delr,r;
2259 
2260  double one = 1.0, two = 2.0;
2261 
2262 
2263 
2264  if(!n)
2265 
2266  return;
2267 
2268 
2269 
2270  for(k = 0; k < n; ++k){
2271 
2272  r = -cos((two*(double)k + one) * dth);
2273 
2274  if(k) r = 0.5*(r + rlast);
2275 
2276 
2277 
2278  for(j = 1; j < STOP; ++j){
2279 
2280  jacobfd(1,&r,&poly, &pder, n, alpha, beta);
2281 
2282 
2283 
2284  for(i = 0, sum = 0.0; i < k; ++i) sum += one/(r - z[i]);
2285 
2286 
2287 
2288  delr = -poly / (pder - sum * poly);
2289 
2290  r += delr;
2291 
2292  if( fabs(delr) < EPS ) break;
2293 
2294  }
2295 
2296  z[k] = r;
2297 
2298  rlast = r;
2299 
2300  }
2301 
2302  return;
2303 
2304  }
STOP
#define STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
Definition: Polylib.cpp:19
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:23
Polylib::jacobfd
void jacobfd(const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
Definition: Polylib.cpp:1920
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 2358 of file Polylib.cpp.

References RecCoeff(), and TriQL().

2360  {
2361 
2362 
2363 
2364  int i,j;
2365 
2366  RecCoeff(n,a,b,alpha,beta);
2367 
2368 
2369 
2370  double **z = new double*[n];
2371 
2372  for(i = 0; i < n; i++)
2373 
2374  {
2375 
2376  z[i] = new double[n];
2377 
2378  for(j = 0; j < n; j++)
2379 
2380  {
2381 
2382  z[i][j] = 0.0;
2383 
2384  }
2385 
2386  }
2387 
2388  for(i = 0; i < n; i++)
2389 
2390  {
2391 
2392  z[i][i] = 1.0;
2393 
2394  }
2395 
2396 
2397 
2398  // find eigenvalues and eigenvectors
2399 
2400  TriQL(n, a, b,z);
2401 
2402 
2403 
2404  return;
2405 
2406  }
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
Definition: Polylib.cpp:2418
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:2534
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 2738 of file Polylib.cpp.

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

2740  {
2741 
2742  int i,j,k,m;
2743 
2744  // Working storage
2745 
2746  int size = (int)floor(n/2.0)+2;
2747 
2748  double *s = new double[size];
2749 
2750  double *t = new double[size];
2751 
2752 
2753 
2754  // Initialize s and t to zero
2755 
2756  for(i = 0; i < size; i++)
2757 
2758  {
2759 
2760  s[i] = 0.0;
2761 
2762  t[i] = 0.0;
2763 
2764  }
2765 
2766 
2767 
2768  t[1] = b[n+1];
2769 
2770  for(m = 0; m <= n-2; m++)
2771 
2772  {
2773 
2774  double u = 0.0;
2775 
2776  for(k = (int)floor((m+1)/2.0); k >= 0; k--)
2777 
2778  {
2779 
2780  int l = m-k;
2781 
2782  u = u+(a[k+n+1]-a[l])*t[k+1] + b[k+n+1]*s[k] - b[l]*s[k+1];
2783 
2784  s[k+1] = u;
2785 
2786  }
2787 
2788 
2789 
2790  // Swap the contents of s and t
2791 
2792  double *hold = s;
2793 
2794  s = t;
2795 
2796  t = hold;
2797 
2798  }
2799 
2800 
2801 
2802 
2803 
2804  for(j = (int)floor(n/2.0); j >= 0; j--)
2805 
2806  {
2807 
2808  s[j+1] = s[j];
2809 
2810  }
2811 
2812 
2813 
2814  for(m = n-1; m <= 2*n-3; m++)
2815 
2816  {
2817 
2818  double u = 0;
2819 
2820  for(k = m+1-n; k <= floor((m-1)/2.0); k++)
2821 
2822  {
2823 
2824  int l = m-k;
2825 
2826  j = n-1-l;
2827 
2828  u = u-(a[k+n+1]-a[l])*t[j+1] - b[k+n+1]*s[j+1] + b[l]*s[j+2];
2829 
2830  s[j+1] = u;
2831 
2832  }
2833 
2834 
2835 
2836  if(m%2 == 0)
2837 
2838  {
2839 
2840  k = m/2;
2841 
2842  a[k+n+1] = a[k] + (s[j+1]-b[k+n+1]*s[j+2])/t[j+2];
2843 
2844 
2845 
2846  }else
2847 
2848  {
2849 
2850  k = (m+1)/2;
2851 
2852  b[k+n+1] = s[j+1]/s[j+2];
2853 
2854  }
2855 
2856 
2857 
2858 
2859 
2860  // Swap the contents of s and t
2861 
2862  double *hold = s;
2863 
2864  s = t;
2865 
2866  t = hold;
2867 
2868  }
2869 
2870 
2871 
2872  a[2*n ] = a[n-1]-b[2*n]*s[1]/t[1];
2873 
2874 
2875 
2876  }
double Polylib::laginterp ( double  z,
int  j,
const double *  zj,
int  np 
)

Definition at line 57 of file Polylib.cpp.

References optdiff().

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

58  {
59  double temp = 1.0;
60  for (int i=0; i<np; i++)
61  {
62  if (j != i)
63  {
64  temp *=optdiff(z,zj[i])/(zj[j]-zj[i]);
65  }
66  }
67  return temp;
68  }
Polylib::optdiff
double optdiff(double xl, double xr)
The following function is used to circumvent/reduce "Subtractive Cancellation" The expression 1/dz is...
Definition: Polylib.cpp:37
double Polylib::optdiff ( double  xl,
double  xr 
)

The following function is used to circumvent/reduce "Subtractive Cancellation" The expression 1/dz is replaced by optinvsub(.,.) Added on 26 April 2017.

Definition at line 37 of file Polylib.cpp.

References L.

Referenced by laginterp().

38  {
39  double m_xln, m_xrn;
40  int m_expn;
41  int m_digits = static_cast<int>(fabs(floor(log10(DBL_EPSILON)))-1);
42 
43  if (fabs(xl-xr)<1.e-4){
44 
45  m_expn = static_cast<int>(floor(log10(fabs(xl-xr))));
46  m_xln = xl*powl(10.0L,-m_expn)-floor(xl*powl(10.0L,-m_expn)); // substract the digits overlap part
47  m_xrn = xr*powl(10.0L,-m_expn)-floor(xl*powl(10.0L,-m_expn)); // substract the common digits overlap part
48  m_xln = round(m_xln*powl(10.0L,m_digits+m_expn)); // git rid of rubbish
49  m_xrn = round(m_xrn*powl(10.0L,m_digits+m_expn));
50 
51  return powl(10.0L,-m_digits)*(m_xln-m_xrn);
52  }else{
53  return (xl-xr);
54  }
55  }
L
NekDouble L
Definition: CompressibleBL.cpp:79
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 2418 of file Polylib.cpp.

References gammaF().

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

2420  {
2421 
2422 
2423 
2424  int i;
2425 
2426  double apb, apbi,a2b2;
2427 
2428 
2429 
2430  if(!n)
2431 
2432  return;
2433 
2434 
2435 
2436  // generate normalised terms
2437 
2438  apb = alpha + beta;
2439 
2440  apbi = 2.0 + apb;
2441 
2442 
2443 
2444  b[0] = pow(2.0,apb+1.0)*gammaF(alpha+1.0)*gammaF(beta+1.0)/gammaF(apbi); //MuZero
2445 
2446  a[0] = (beta-alpha)/apbi;
2447 
2448  b[1] = (4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));
2449 
2450 
2451 
2452  a2b2 = beta*beta-alpha*alpha;
2453 
2454 
2455 
2456  for(i = 1; i < n-1; i++){
2457 
2458  apbi = 2.0*(i+1) + apb;
2459 
2460  a[i] = a2b2/((apbi-2.0)*apbi);
2461 
2462  b[i+1] = (4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
2463 
2464  ((apbi*apbi-1)*apbi*apbi));
2465 
2466  }
2467 
2468 
2469 
2470  apbi = 2.0*n + apb;
2471 
2472  a[n-1] = a2b2/((apbi-2.0)*apbi);
2473 
2474 
2475 
2476  }
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2174
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 2534 of file Polylib.cpp.

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

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

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

Gauss-Jacobi zeros and weights.

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

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

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

Definition at line 142 of file Polylib.cpp.

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

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

146  {
147 
148  register int i;
149 
150  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
151 
152 
153 
154  jacobz (np,z,alpha,beta);
155 
156  jacobd (np,z,w,np,alpha,beta);
157 
158 
159 
160  fac = pow(two,apb + one)*gammaF(alpha + np + one)*gammaF(beta + np + one);
161 
162  fac /= gammaF(np + one)*gammaF(apb + np + one);
163 
164 
165 
166  for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]*(one-z[i]*z[i]));
167 
168 
169 
170  return;
171 
172  }
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:2120
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:80
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2174
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 420 of file Polylib.cpp.

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

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

424  {
425 
426 
427 
428  int np = (npt-1)/2;
429 
430 
431 
432  int i,j;
433 
434 
435 
436  // number of kronrod points associated with the np gauss rule
437 
438  int kpoints = 2*np + 1;
439 
440 
441 
442  // Define the number of required recurrence coefficents
443 
444  int ncoeffs = (int)floor(3.0*(np+1)/2);
445 
446 
447 
448  // Define arrays for the recurrence coefficients
449 
450  // We will use these arrays for the Kronrod results too, hence the
451 
452  // reason for the size of the arrays
453 
454  double *a = new double[kpoints];
455 
456  double *b = new double[kpoints];
457 
458 
459 
460  // Initialize a and b to zero
461 
462  for(i = 0; i < kpoints; i++)
463 
464  {
465 
466  a[i] = 0.0;
467 
468  b[i] = 0.0;
469 
470  }
471 
472 
473 
474  // Call the routine to calculate the recurrence coefficients
475 
476  RecCoeff(ncoeffs,a,b,alpha,beta);
477 
478 
479 
480  // Call the routine to caluclate the jacobi-Kronrod matrix
481 
482  JKMatrix(np,a,b);
483 
484 
485 
486  // Set up the identity matrix
487 
488  double** zmatrix = new double*[kpoints];
489 
490  for(i = 0; i < kpoints; i++)
491 
492  {
493 
494  zmatrix[i] = new double[kpoints];
495 
496  for(j = 0; j < kpoints; j++)
497 
498  {
499 
500  zmatrix[i][j] = 0.0;
501 
502  }
503 
504  }
505 
506  for(i = 0; i < kpoints; i++)
507 
508  {
509 
510  zmatrix[i][i] = 1.0;
511 
512  }
513 
514 
515 
516  // Calculte the points and weights
517 
518  TriQL(kpoints, a, b, zmatrix);
519 
520 
521 
522  for(i = 0; i < kpoints; i++)
523 
524  {
525 
526  z[i] = a[i];
527 
528  w[i] = b[i];
529 
530  }
531 
532 
533 
534 
535 
536  }
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
Definition: Polylib.cpp:2418
Polylib::JKMatrix
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the.
Definition: Polylib.cpp:2738
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:2534
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 346 of file Polylib.cpp.

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

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

350  {
351 
352 
353 
354  if( np == 1 ){
355 
356  z[0] = 0.0;
357 
358  w[0] = 2.0;
359 
360  }
361 
362  else{
363 
364  register int i;
365 
366  double fac, one = 1.0, apb = alpha + beta, two = 2.0;
367 
368 
369 
370  z[0] = -one;
371 
372  z[np-1] = one;
373 
374  jacobz (np-2,z + 1,alpha + one,beta + one);
375 
376  jacobfd (np,z,w,NULL,np-1,alpha,beta);
377 
378 
379 
380  fac = pow(two,apb + 1)*gammaF(alpha + np)*gammaF(beta + np);
381 
382  fac /= (np-1)*gammaF(np)*gammaF(alpha + beta + np + one);
383 
384 
385 
386  for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]);
387 
388  w[0] *= (beta + one);
389 
390  w[np-1] *= (alpha + one);
391 
392  }
393 
394 
395 
396  return;
397 
398  }
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:80
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:1920
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2174
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 198 of file Polylib.cpp.

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

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

202  {
203 
204 
205 
206  if(np == 1){
207 
208  z[0] = 0.0;
209 
210  w[0] = 2.0;
211 
212  }
213 
214  else{
215 
216  register int i;
217 
218  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
219 
220 
221 
222  z[0] = -one;
223 
224  jacobz (np-1,z+1,alpha,beta+1);
225 
226  jacobfd (np,z,w,NULL,np-1,alpha,beta);
227 
228 
229 
230  fac = pow(two,apb)*gammaF(alpha + np)*gammaF(beta + np);
231 
232  fac /= gammaF(np)*(beta + np)*gammaF(apb + np + 1);
233 
234 
235 
236  for(i = 0; i < np; ++i) w[i] = fac*(1-z[i])/(w[i]*w[i]);
237 
238  w[0] *= (beta + one);
239 
240  }
241 
242 
243 
244  return;
245 
246  }
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:80
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:1920
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2174
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 274 of file Polylib.cpp.

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

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

278  {
279 
280 
281 
282  if(np == 1){
283 
284  z[0] = 0.0;
285 
286  w[0] = 2.0;
287 
288  }
289 
290  else{
291 
292  register int i;
293 
294  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
295 
296 
297 
298  jacobz (np-1,z,alpha+1,beta);
299 
300  z[np-1] = one;
301 
302  jacobfd (np,z,w,NULL,np-1,alpha,beta);
303 
304 
305 
306  fac = pow(two,apb)*gammaF(alpha + np)*gammaF(beta + np);
307 
308  fac /= gammaF(np)*(alpha + np)*gammaF(apb + np + 1);
309 
310 
311 
312  for(i = 0; i < np; ++i) w[i] = fac*(1+z[i])/(w[i]*w[i]);
313 
314  w[np-1] *= (alpha + one);
315 
316  }
317 
318 
319 
320  return;
321 
322  }
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:80
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:1920
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2174
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 746 of file Polylib.cpp.

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

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

750  {
751 
752 
753 
754  int np = (npt+1)/2;
755 
756 
757 
758  if(np < 4)
759 
760  {
761 
762  fprintf (stderr,"too few points in formula\n");
763 
764  return;
765 
766  }
767 
768 
769 
770  double endl = -1;
771 
772  double endr = 1;
773 
774  int i,j;
775 
776 
777 
778  // number of kronrod points associated with the np gauss rule
779 
780  int kpoints = 2*np-1;
781 
782 
783 
784  // Define the number of required recurrence coefficents
785 
786  int ncoeffs = (int)ceil(3.0*np/2)-1;
787 
788 
789 
790  // Define arrays for the recurrence coefficients
791 
792  double *a = new double[ncoeffs+1];
793 
794  double *b = new double[ncoeffs+1];
795 
796 
797 
798  // Initialize a and b to zero
799 
800  for(i = 0; i < ncoeffs+1; i++)
801 
802  {
803 
804  a[i] = 0.0;
805 
806  b[i] = 0.0;
807 
808  }
809 
810 
811 
812  // Call the routine to calculate the recurrence coefficients
813 
814  RecCoeff(ncoeffs,a,b,alpha,beta);
815 
816 
817 
818 
819 
820  double* a0 = new double[ncoeffs];
821 
822  double* b0 = new double[ncoeffs];
823 
824 
825 
826  chri1(ncoeffs,a,b,a0,b0,endl);
827 
828 
829 
830  double* a1 = new double[ncoeffs-1];
831 
832  double* b1 = new double[ncoeffs-1];
833 
834 
835 
836  chri1(ncoeffs-1,a0,b0,a1,b1,endr);
837 
838 
839 
840 
841 
842  double s = b1[0]/fabs(b1[0]);
843 
844  b1[0] = s*b1[0];
845 
846 
847 
848  // Finding the 2*np-1 gauss-kronrod points
849 
850  double* z1 = new double[2*np-3];
851 
852  double* w1 = new double[2*np-3];
853 
854  for(i = 0; i < ncoeffs; i++)
855 
856  {
857 
858  z1[i] = a1[i];
859 
860  w1[i] = b1[i];
861 
862  }
863 
864  JKMatrix(np-2,z1,w1);
865 
866  // Set up the identity matrix
867 
868  double** zmatrix = new double*[2*np-3];
869 
870  for(i = 0; i < 2*np-3; i++)
871 
872  {
873 
874  zmatrix[i] = new double[2*np-3];
875 
876  for(j = 0; j < 2*np-3; j++)
877 
878  {
879 
880  zmatrix[i][j] = 0.0;
881 
882  }
883 
884  }
885 
886  for(i = 0; i < 2*np-3; i++)
887 
888  {
889 
890  zmatrix[i][i] = 1.0;
891 
892  }
893 
894 
895 
896  // Calculate the points and weights
897 
898  TriQL(2*np-3, z1, w1, zmatrix);
899 
900 
901 
902  double sumW1 = 0.0;
903 
904  double sumW1Z1 = 0.0;
905 
906  for(i = 0; i < 2*np-3; i++)
907 
908  {
909 
910  w1[i] = s*w1[i]/(z1[i]-endl)/(z1[i]-endr);
911 
912  sumW1 += w1[i];
913 
914  sumW1Z1 += z1[i]*w1[i];
915 
916  }
917 
918 
919 
920  double c0 = b[0]-sumW1;
921 
922  double c1 = a[0]*b[0]-sumW1Z1;
923 
924 
925 
926  z[0] = endl;
927 
928  z[2*np-2] = endr;
929 
930  w[0] = (c0*endr-c1)/(endr-endl);
931 
932  w[2*np-2] = (c1-c0*endl)/(endr-endl);
933 
934 
935 
936  for(i = 1; i < kpoints-1; i++)
937 
938  {
939 
940  z[i] = z1[i-1];
941 
942  w[i] = w1[i-1];
943 
944  }
945 
946  }
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
Definition: Polylib.cpp:2418
Polylib::chri1
void chri1(int, double *, double *, double *, double *, double)
Given a weight function through the first n+1.
Definition: Polylib.cpp:2900
Polylib::JKMatrix
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the.
Definition: Polylib.cpp:2738
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:2534
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 554 of file Polylib.cpp.

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

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

558  {
559 
560 
561 
562  int np = npt/2;
563 
564 
565 
566  if(np < 2)
567 
568  {
569 
570  fprintf(stderr,"too few points in formula\n");
571 
572  return;
573 
574  }
575 
576 
577 
578  double end0 = -1;
579 
580 
581 
582  int i,j;
583 
584 
585 
586  // number of kronrod points associated with the np gauss rule
587 
588  int kpoints = 2*np;
589 
590 
591 
592  // Define the number of required recurrence coefficents
593 
594  int ncoeffs = (int)ceil(3.0*np/2);
595 
596 
597 
598  // Define arrays for the recurrence coefficients
599 
600  double *a = new double[ncoeffs+1];
601 
602  double *b = new double[ncoeffs+1];
603 
604 
605 
606  // Initialize a and b to zero
607 
608  for(i = 0; i < ncoeffs+1; i++)
609 
610  {
611 
612  a[i] = 0.0;
613 
614  b[i] = 0.0;
615 
616  }
617 
618 
619 
620  // Call the routine to calculate the recurrence coefficients
621 
622  RecCoeff(ncoeffs,a,b,alpha,beta);
623 
624 
625 
626  double* a0 = new double[ncoeffs];
627 
628  double* b0 = new double[ncoeffs];
629 
630 
631 
632  chri1(ncoeffs,a,b,a0,b0,end0);
633 
634 
635 
636  double s = b0[0]/fabs(b0[0]);
637 
638  b0[0] = s*b0[0];
639 
640 
641 
642  // Finding the 2*np-1 gauss-kronrod points
643 
644  double* z1 = new double[2*np-1];
645 
646  double* w1 = new double[2*np-1];
647 
648  for(i = 0; i < ncoeffs; i++)
649 
650  {
651 
652  z1[i] = a0[i];
653 
654  w1[i] = b0[i];
655 
656  }
657 
658  JKMatrix(np-1,z1,w1);
659 
660  // Set up the identity matrix
661 
662  double** zmatrix = new double*[2*np-1];
663 
664  for(i = 0; i < 2*np-1; i++)
665 
666  {
667 
668  zmatrix[i] = new double[2*np-1];
669 
670  for(j = 0; j < 2*np-1; j++)
671 
672  {
673 
674  zmatrix[i][j] = 0.0;
675 
676  }
677 
678  }
679 
680  for(i = 0; i < 2*np-1; i++)
681 
682  {
683 
684  zmatrix[i][i] = 1.0;
685 
686  }
687 
688 
689 
690  // Calculate the points and weights
691 
692  TriQL(2*np-1, z1, w1, zmatrix);
693 
694 
695 
696  double sumW1 = 0.0;
697 
698  for(i = 0; i < 2*np-1; i++)
699 
700  {
701 
702  w1[i] = s*w1[i]/(z1[i]-end0);
703 
704  sumW1 += w1[i];
705 
706  }
707 
708 
709 
710  z[0] = end0;
711 
712  w[0] = b[0]- sumW1;
713 
714  for(i = 1; i < kpoints; i++)
715 
716  {
717 
718  z[i] = z1[i-1];
719 
720  w[i] = w1[i-1];
721 
722  }
723 
724 
725 
726 
727 
728  }
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
Definition: Polylib.cpp:2418
Polylib::chri1
void chri1(int, double *, double *, double *, double *, double)
Given a weight function through the first n+1.
Definition: Polylib.cpp:2900
Polylib::JKMatrix
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the.
Definition: Polylib.cpp:2738
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:2534
Generated on Sat Nov 25 2017 16:19:13 for Nektar++ by   doxygen 1.8.8