Nektar++
Functions
Polylib Namespace Reference

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

Functions

double optdiff (double xl, double xr)
 The following function is used to circumvent/reduce "Subtractive Cancellation" The expression 1/dz is replaced by optinvsub(.,.) Added on 26 April 2017. More...
 
double laginterp (double z, int j, const double *zj, int np)
 
double laginterpderiv (double z, int k, const double *zj, int np)
 
static void Jacobz (const int n, double *z, const double alpha, const double beta)
 Calculate the n zeros, z, of the Jacobi polynomial, i.e. \( P_n^{\alpha,\beta}(z) = 0 \). More...
 
static void RecCoeff (const int n, double *a, double *b, const double alpha, const double beta)
 The routine finds the recurrence coefficients a and b of the orthogonal polynomials. More...
 
static void TriQL (const int n, double *d, double *e, double **z)
 QL algorithm for symmetric tridiagonal matrix. More...
 
void JKMatrix (int n, double *a, double *b)
 Calcualtes the Jacobi-kronrod matrix by determining the a and coefficients. More...
 
void chri1 (int n, double *a, double *b, double *a0, double *b0, double z)
 Given a weight function \(w(t)\) through the first n+1 coefficients a and b of its orthogonal polynomials this routine generates the first n recurrence coefficients for the orthogonal polynomials relative to the modified weight function \((t-z)w(t)\). More...
 
void zwgj (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Jacobi zeros and weights. More...
 
void zwgrjm (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Radau-Jacobi zeros and weights with end point at z=-1. More...
 
void zwgrjp (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Radau-Jacobi zeros and weights with end point at z=1. More...
 
void zwglj (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1. More...
 
void zwgk (double *z, double *w, const int npt, const double alpha, const double beta)
 Gauss-Kronrod-Jacobi zeros and weights. More...
 
void zwrk (double *z, double *w, const int npt, const double alpha, const double beta)
 Gauss-Radau-Kronrod-Jacobi zeros and weights. More...
 
void zwlk (double *z, double *w, const int npt, const double alpha, const double beta)
 Gauss-Lobatto-Kronrod-Jacobi zeros and weights. More...
 
void Qg (double *Q, const double *z, const int np)
 Compute the Integration Matrix. More...
 
void Dgj (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix and its transpose associated with the Gauss-Jacobi zeros. More...
 
void Dgrjm (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix and its transpose associated with the Gauss-Radau-Jacobi zeros with a zero at z=-1. More...
 
void Dgrjp (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix associated with the Gauss-Radau-Jacobi zeros with a zero at z=1. More...
 
void Dglj (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix associated with the Gauss-Lobatto-Jacobi zeros. More...
 
double hgj (const int i, const double z, const double *zgj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through the np Gauss-Jacobi points zgj at the arbitrary location z. More...
 
double hgrjm (const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at the arbitrary location z. This routine assumes zgrj includes the point -1. More...
 
double hgrjp (const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at the arbitrary location z. This routine assumes zgrj includes the point +1. More...
 
double hglj (const int i, const double z, const double *zglj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through the np Gauss-Lobatto-Jacobi points zgrj at the arbitrary location z. More...
 
void Imgj (double *im, const double *zgj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Jacobi points to an arbitrary distribution at points zm. More...
 
void Imgrjm (double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Radau-Jacobi points (including z=-1) to an arbitrary distrubtion at points zm. More...
 
void Imgrjp (double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Radau-Jacobi points (including z=1) to an arbitrary distrubtion at points zm. More...
 
void Imglj (double *im, const double *zglj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Lobatto-Jacobi points to an arbitrary distrubtion at points zm. More...
 
void polycoeffs (double *c, const double *z, const int i, const int np)
 Compute the coefficients of Lagrange interpolation polynomials. 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...
 
double gammaF (const double x)
 Calculate the Gamma function , \( \Gamma(n)\), for integer values and halves. More...
 
double gammaFracGammaF (const int x, const double alpha, const int y, const double beta)
 Calculate fraction of two Gamma functions, \( \Gamma(x+\alpha)/\Gamma(y+\beta) \), for integer values and halves. More...
 
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...
 
void JacZeros (const int n, double *a, double *b, const double alpha, const double beta)
 Zero and Weight determination through the eigenvalues and eigenvectors of a tridiagonal matrix from the three term recurrence relationship. More...
 

Detailed Description

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

Function Documentation

◆ chri1()

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

Given a weight function \(w(t)\) through the first n+1 coefficients a and b of its orthogonal polynomials this routine generates the first n recurrence coefficients for the orthogonal polynomials relative to the modified weight function \((t-z)w(t)\).

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

Definition at line 1956 of file Polylib.cpp.

1957{
1958
1959 double q = ceil(3.0 * n / 2);
1960 int size = (int)q + 1;
1961 if (size < n + 1)
1962 {
1963 fprintf(stderr, "input arrays a and b are too short\n");
1964 }
1965 double *r = new double[n + 1];
1966 r[0] = z - a[0];
1967 r[1] = z - a[1] - b[1] / r[0];
1968 a0[0] = a[1] + r[1] - r[0];
1969 b0[0] = -r[0] * b[0];
1970
1971 if (n == 1)
1972 {
1973 delete[] r;
1974 return;
1975 }
1976 int k = 0;
1977 for (k = 1; k < n; k++)
1978 {
1979 r[k + 1] = z - a[k + 1] - b[k + 1] / r[k];
1980 a0[k] = a[k + 1] + r[k + 1] - r[k];
1981 b0[k] = b[k] * r[k] / r[k - 1];
1982 }
1983 delete[] r;
1984}
Nektar::UnitTests::z
std::vector< double > z(NPUPPER)
Nektar::UnitTests::q
std::vector< double > q(NPUPPER *NPUPPER)

References Nektar::UnitTests::q(), and Nektar::UnitTests::z().

Referenced by zwlk(), and zwrk().

◆ Dgj()

void Polylib::Dgj ( double *  D,
const double *  z,
const int  np,
const double  alpha,
const double  beta 
)

Compute the Derivative Matrix and its transpose associated with the Gauss-Jacobi zeros.

  • Compute the derivative matrix, d, and its transpose, dt, associated with the n_th order Lagrangian interpolants through the np Gauss-Jacobi points z such that
    \( \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \)

Definition at line 653 of file Polylib.cpp.

655{
656
657 double one = 1.0, two = 2.0;
658
659 if (np <= 0)
660 {
661 D[0] = 0.0;
662 }
663 else
664 {
665 int i, j;
666 double *pd;
667
668 pd = (double *)malloc(np * sizeof(double));
669 jacobd(np, z, pd, np, alpha, beta);
670
671 for (i = 0; i < np; i++)
672 {
673 for (j = 0; j < np; j++)
674 {
675
676 if (i != j)
677 {
678 D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
679 }
680 else
681 {
682 D[i * np + j] =
683 (alpha - beta + (alpha + beta + two) * z[j]) /
684 (two * (one - z[j] * z[j]));
685 }
686 }
687 }
688 free(pd);
689 }
690 return;
691}
Nektar::LibUtilities::beta
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:59
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:1378

References Nektar::LibUtilities::beta, jacobd(), and Nektar::UnitTests::z().

Referenced by Nektar::UnitTests::BOOST_AUTO_TEST_CASE(), and Nektar::LibUtilities::GaussPoints::v_CalculateDerivMatrix().

◆ Dglj()

void Polylib::Dglj ( double *  D,
const double *  z,
const int  np,
const double  alpha,
const double  beta 
)

Compute the Derivative Matrix associated with the Gauss-Lobatto-Jacobi zeros.

  • Compute the derivative matrix, d, associated with the n_th order Lagrange interpolants through the np Gauss-Lobatto-Jacobi points z such that
    \( \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \)

Definition at line 833 of file Polylib.cpp.

835{
836 if (np <= 1)
837 {
838 D[0] = 0.0;
839 }
840 else
841 {
842 int i, j;
843 double one = 1.0, two = 2.0;
844 double *pd;
845
846 pd = (double *)malloc(np * sizeof(double));
847
848 pd[0] = two * pow(-one, np) * gammaFracGammaF(np, beta, np - 1, 0.0);
849 pd[0] /= gammaF(beta + two);
850 jacobd(np - 2, z + 1, pd + 1, np - 2, alpha + 1, beta + 1);
851 for (i = 1; i < np - 1; ++i)
852 {
853 pd[i] *= (one - z[i] * z[i]);
854 }
855 pd[np - 1] = -two * gammaFracGammaF(np, alpha, np - 1, 0.0);
856 pd[np - 1] /= gammaF(alpha + two);
857
858 for (i = 0; i < np; i++)
859 {
860 for (j = 0; j < np; j++)
861 {
862 if (i != j)
863 {
864 D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
865 }
866 else
867 {
868 if (j == 0)
869 {
870 D[i * np + j] =
871 (alpha - (np - 1) * (np + alpha + beta)) /
872 (two * (beta + two));
873 }
874 else if (j == np - 1)
875 {
876 D[i * np + j] =
877 -(beta - (np - 1) * (np + alpha + beta)) /
878 (two * (alpha + two));
879 }
880 else
881 {
882 D[i * np + j] = (alpha - beta + (alpha + beta) * z[j]) /
883 (two * (one - z[j] * z[j]));
884 }
885 }
886 }
887 }
888 free(pd);
889 }
890
891 return;
892}
Polylib::gammaF
double gammaF(const double x)
Calculate the Gamma function , , for integer values and halves.
Definition: Polylib.cpp:1413
Polylib::gammaFracGammaF
double gammaFracGammaF(const int x, const double alpha, const int y, const double beta)
Calculate fraction of two Gamma functions, , for integer values and halves.
Definition: Polylib.cpp:1468

References Nektar::LibUtilities::beta, gammaF(), gammaFracGammaF(), jacobd(), and Nektar::UnitTests::z().

Referenced by Nektar::UnitTests::BOOST_AUTO_TEST_CASE(), and Nektar::LibUtilities::GaussPoints::v_CalculateDerivMatrix().

◆ Dgrjm()

void Polylib::Dgrjm ( double *  D,
const double *  z,
const int  np,
const double  alpha,
const double  beta 
)

Compute the Derivative Matrix and its transpose associated with the Gauss-Radau-Jacobi zeros with a zero at z=-1.

  • Compute the derivative matrix, d, associated with the n_th order Lagrangian interpolants through the np Gauss-Radau-Jacobi points z such that
    \( \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \)

Definition at line 704 of file Polylib.cpp.

706{
707
708 if (np <= 0)
709 {
710 D[0] = 0.0;
711 }
712 else
713 {
714 int i, j;
715 double one = 1.0, two = 2.0;
716 double *pd;
717
718 pd = (double *)malloc(np * sizeof(double));
719
720 pd[0] = pow(-one, np - 1) * gammaFracGammaF(np + 1, beta, np, 0.0);
721 pd[0] /= gammaF(beta + two);
722 jacobd(np - 1, z + 1, pd + 1, np - 1, alpha, beta + 1);
723 for (i = 1; i < np; ++i)
724 {
725 pd[i] *= (1 + z[i]);
726 }
727
728 for (i = 0; i < np; i++)
729 {
730 for (j = 0; j < np; j++)
731 {
732 if (i != j)
733 {
734 D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
735 }
736 else
737 {
738 if (j == 0)
739 {
740 D[i * np + j] = -(np + alpha + beta + one) *
741 (np - one) / (two * (beta + two));
742 }
743 else
744 {
745 D[i * np + j] =
746 (alpha - beta + one + (alpha + beta + one) * z[j]) /
747 (two * (one - z[j] * z[j]));
748 }
749 }
750 }
751 }
752 free(pd);
753 }
754
755 return;
756}

References Nektar::LibUtilities::beta, gammaF(), gammaFracGammaF(), jacobd(), and Nektar::UnitTests::z().

Referenced by Nektar::UnitTests::BOOST_AUTO_TEST_CASE(), and Nektar::LibUtilities::GaussPoints::v_CalculateDerivMatrix().

◆ Dgrjp()

void Polylib::Dgrjp ( double *  D,
const double *  z,
const int  np,
const double  alpha,
const double  beta 
)

Compute the Derivative Matrix associated with the Gauss-Radau-Jacobi zeros with a zero at z=1.

  • Compute the derivative matrix, d, associated with the n_th order Lagrangian interpolants through the np Gauss-Radau-Jacobi points z such that
    \( \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \)

Definition at line 768 of file Polylib.cpp.

770{
771
772 if (np <= 0)
773 {
774 D[0] = 0.0;
775 }
776 else
777 {
778 int i, j;
779 double one = 1.0, two = 2.0;
780 double *pd;
781
782 pd = (double *)malloc(np * sizeof(double));
783
784 jacobd(np - 1, z, pd, np - 1, alpha + 1, beta);
785 for (i = 0; i < np - 1; ++i)
786 {
787 pd[i] *= (1 - z[i]);
788 }
789 pd[np - 1] = -gammaFracGammaF(np + 1, alpha, np, 0.0);
790 pd[np - 1] /= gammaF(alpha + two);
791
792 for (i = 0; i < np; i++)
793 {
794 for (j = 0; j < np; j++)
795 {
796 if (i != j)
797 {
798 D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
799 }
800 else
801 {
802 if (j == np - 1)
803 {
804 D[i * np + j] = (np + alpha + beta + one) * (np - one) /
805 (two * (alpha + two));
806 }
807 else
808 {
809 D[i * np + j] =
810 (alpha - beta - one + (alpha + beta + one) * z[j]) /
811 (two * (one - z[j] * z[j]));
812 }
813 }
814 }
815 }
816 free(pd);
817 }
818
819 return;
820}

References Nektar::LibUtilities::beta, gammaF(), gammaFracGammaF(), jacobd(), and Nektar::UnitTests::z().

Referenced by Nektar::UnitTests::BOOST_AUTO_TEST_CASE(), and Nektar::LibUtilities::GaussPoints::v_CalculateDerivMatrix().

◆ gammaF()

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

1414{
1415 double gamma = 1.0;
1416
1417 if (x == -0.5)
1418 {
1419 gamma = -2.0 * sqrt(M_PI);
1420 }
1421 else if (!x)
1422 {
1423 return gamma;
1424 }
1425 else if ((x - (int)x) == 0.5)
1426 {
1427 int n = (int)x;
1428 double tmp = x;
1429
1430 gamma = sqrt(M_PI);
1431 while (n--)
1432 {
1433 tmp -= 1.0;
1434 gamma *= tmp;
1435 }
1436 }
1437 else if ((x - (int)x) == 0.0)
1438 {
1439 int n = (int)x;
1440 double tmp = x;
1441
1442 while (--n)
1443 {
1444 tmp -= 1.0;
1445 gamma *= tmp;
1446 }
1447 }
1448 else
1449 {
1450 fprintf(stderr, "%lf is not of integer or half order\n", x);
1451 }
1452 return gamma;
1453}
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294

References tinysimd::sqrt().

Referenced by Nektar::UnitTests::BOOST_AUTO_TEST_CASE(), Dglj(), Dgrjm(), Dgrjp(), and RecCoeff().

◆ gammaFracGammaF()

double Polylib::gammaFracGammaF ( const int  x,
const double  alpha,
const int  y,
const double  beta 
)

Calculate fraction of two Gamma functions, \( \Gamma(x+\alpha)/\Gamma(y+\beta) \), for integer values and halves.

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

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

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

Attempts simplification in cases like \( \Gamma(x+1)/\Gamma(x-1) = x*(x-1) \)

Definition at line 1468 of file Polylib.cpp.

1470{
1471 double gamma = 1.0;
1472 double halfa = fabs(alpha - int(alpha));
1473 double halfb = fabs(beta - int(beta));
1474 if (halfa == 0.0 && halfb == 0.0) // integer value
1475 {
1476 int X = x + alpha;
1477 int Y = y + beta;
1478 if (X > Y)
1479 {
1480 for (int tmp = X - 1; tmp > Y - 1; tmp -= 1)
1481 {
1482 gamma *= tmp;
1483 }
1484 }
1485 else if (Y > X)
1486 {
1487 for (int tmp = Y - 1; tmp > X - 1; tmp -= 1)
1488 {
1489 gamma *= tmp;
1490 }
1491 gamma = 1. / gamma;
1492 }
1493 }
1494 else if (halfa == 0.5 && halfb == 0.5) // both are halves
1495 {
1496 double X = x + alpha;
1497 double Y = y + beta;
1498 if (X > Y)
1499 {
1500 for (int tmp = int(X); tmp > int(Y); tmp -= 1)
1501 {
1502 gamma *= tmp - 0.5;
1503 }
1504 }
1505 else if (Y > X)
1506 {
1507 for (int tmp = int(Y); tmp > int(X); tmp -= 1)
1508 {
1509 gamma *= tmp - 0.5;
1510 }
1511 gamma = 1. / gamma;
1512 }
1513 }
1514 else
1515 {
1516 double X = x + alpha;
1517 double Y = y + beta;
1518 while (X > 1 || Y > 1)
1519 {
1520 if (X > 1)
1521 {
1522 gamma *= X - 1.;
1523 X -= 1.;
1524 }
1525 if (Y > 1)
1526 {
1527 gamma /= Y - 1.;
1528 Y -= 1.;
1529 }
1530 }
1531 if (X == 0.5)
1532 {
1533 gamma *= sqrt(M_PI);
1534 }
1535 else if (Y == 0.5)
1536 {
1537 gamma /= sqrt(M_PI);
1538 }
1539 else
1540 {
1541 fprintf(stderr, "%lf or %lf is not of integer or half order\n", X,
1542 Y);
1543 }
1544 }
1545
1546 return gamma;
1547}

References Nektar::LibUtilities::beta, and tinysimd::sqrt().

Referenced by Nektar::UnitTests::BOOST_AUTO_TEST_CASE(), Dglj(), Dgrjm(), Dgrjp(), zwgj(), zwglj(), zwgrjm(), and zwgrjp().

◆ hgj()

double Polylib::hgj ( const int  i,
const double  z,
const double *  zgj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through the np Gauss-Jacobi points zgj at the arbitrary location z.

  • \( -1 \leq z \leq 1 \)
  • Uses the defintion of the Lagrangian interpolant:
    % \( \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{P_{np}^{\alpha,\beta}(z)} {[P_{np}^{\alpha,\beta}(z_j)]^\prime (z-z_j)} & \mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} \)

Definition at line 914 of file Polylib.cpp.

917{
918 double zi, dz;
919
920 zi = *(zgj + i);
921 dz = z - zi;
922 if (fabs(dz) < EPS)
923 {
924 return 1.0;
925 }
926
927 return laginterp(z, i, zgj, np);
928}
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:45
Polylib::laginterp
double laginterp(double z, int j, const double *zj, int np)
Definition: Polylib.cpp:85

References EPS, laginterp(), and Nektar::UnitTests::z().

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

◆ hglj()

double Polylib::hglj ( const int  i,
const double  z,
const double *  zglj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through the np Gauss-Lobatto-Jacobi points zgrj at the arbitrary location z.

  • \( -1 \leq z \leq 1 \)
  • Uses the defintion of the Lagrangian interpolant:
    % \( \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{(1-z^2) P_{np-2}^{\alpha+1,\beta+1}(z)} {((1-z^2_j) [P_{np-2}^{\alpha+1,\beta+1}(z_j)]^\prime - 2 z_j P_{np-2}^{\alpha+1,\beta+1}(z_j) ) (z-z_j)}&\mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} \)

Definition at line 1025 of file Polylib.cpp.

1028{
1029 double zi, dz;
1030
1031 zi = *(zglj + i);
1032 dz = z - zi;
1033 if (fabs(dz) < EPS)
1034 {
1035 return 1.0;
1036 }
1037
1038 return laginterp(z, i, zglj, np);
1039}

References EPS, laginterp(), and Nektar::UnitTests::z().

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

◆ hgrjm()

double Polylib::hgrjm ( const int  i,
const double  z,
const double *  zgrj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at the arbitrary location z. This routine assumes zgrj includes the point -1.

  • \( -1 \leq z \leq 1 \)
  • Uses the defintion of the Lagrangian interpolant:
    % \( \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{(1+z) P_{np-1}^{\alpha,\beta+1}(z)} {((1+z_j) [P_{np-1}^{\alpha,\beta+1}(z_j)]^\prime + P_{np-1}^{\alpha,\beta+1}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} \)

Definition at line 951 of file Polylib.cpp.

954{
955 double zi, dz;
956
957 zi = *(zgrj + i);
958 dz = z - zi;
959 if (fabs(dz) < EPS)
960 {
961 return 1.0;
962 }
963
964 return laginterp(z, i, zgrj, np);
965}

References EPS, laginterp(), and Nektar::UnitTests::z().

Referenced by Imgrjm().

◆ hgrjp()

double Polylib::hgrjp ( const int  i,
const double  z,
const double *  zgrj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at the arbitrary location z. This routine assumes zgrj includes the point +1.

  • \( -1 \leq z \leq 1 \)
  • Uses the defintion of the Lagrangian interpolant:
    % \( \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{(1-z) P_{np-1}^{\alpha+1,\beta}(z)} {((1-z_j) [P_{np-1}^{\alpha+1,\beta}(z_j)]^\prime - P_{np-1}^{\alpha+1,\beta}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} \)

Definition at line 988 of file Polylib.cpp.

991{
992 double zi, dz;
993
994 zi = *(zgrj + i);
995 dz = z - zi;
996 if (fabs(dz) < EPS)
997 {
998 return 1.0;
999 }
1000
1001 return laginterp(z, i, zgrj, np);
1002}

References EPS, laginterp(), and Nektar::UnitTests::z().

Referenced by Imgrjp().

◆ ImagBesselComp()

std::complex< Nektar::NekDouble > Polylib::ImagBesselComp ( int  n,
std::complex< Nektar::NekDouble >  y 
)

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

Returns a complex double

Definition at line 1559 of file Polylib.cpp.

1561{
1562 std::complex<Nektar::NekDouble> z(1.0, 0.0);
1563 std::complex<Nektar::NekDouble> zbes(1.0, 0.0);
1564 std::complex<Nektar::NekDouble> zarg;
1565 Nektar::NekDouble tol = 1e-15;
1566 int maxit = 10000;
1567 int i = 1;
1568
1569 zarg = -0.25 * y * y;
1570
1571 while (abs(z) > tol && i <= maxit)
1572 {
1573 z = z * (1.0 / i / (i + n) * zarg);
1574 if (abs(z) <= tol)
1575 {
1576 break;
1577 }
1578 zbes = zbes + z;
1579 i++;
1580 }
1581 zarg = 0.5 * y;
1582 for (i = 1; i <= n; i++)
1583 {
1584 zbes = zbes * zarg;
1585 }
1586 return zbes;
1587}
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:298

References tinysimd::abs(), and Nektar::UnitTests::z().

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

◆ Imgj()

void Polylib::Imgj ( double *  im,
const double *  zgj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Jacobi points to an arbitrary distribution at points zm.

  • Computes the one-dimensional interpolation matrix, im, to interpolate a function from at Gauss-Jacobi distribution of nz zeros zgrj to an arbitrary distribution of mz points zm, i.e.
    \( u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \)

Definition at line 1054 of file Polylib.cpp.

1056{
1057 double zp;
1058 int i, j;
1059
1060 for (i = 0; i < nz; ++i)
1061 {
1062 for (j = 0; j < mz; ++j)
1063 {
1064 zp = zm[j];
1065 im[i * mz + j] = hgj(i, zp, zgj, nz, alpha, beta);
1066 }
1067 }
1068
1069 return;
1070}
Polylib::hgj
double hgj(const int i, const double z, const double *zgj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Jacobi points zgj at the ar...
Definition: Polylib.cpp:914

References Nektar::LibUtilities::beta, and hgj().

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

◆ Imglj()

void Polylib::Imglj ( double *  im,
const double *  zglj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Lobatto-Jacobi points to an arbitrary distrubtion at points zm.

  • Computes the one-dimensional interpolation matrix, im, to interpolate a function from at Gauss-Lobatto-Jacobi distribution of nz zeros zgrj (where zgrj[0]=-1) to an arbitrary distribution of mz points zm, i.e.
    \( u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \)

Definition at line 1147 of file Polylib.cpp.

1149{
1150 double zp;
1151 int i, j;
1152
1153 for (i = 0; i < nz; i++)
1154 {
1155 for (j = 0; j < mz; j++)
1156 {
1157 zp = zm[j];
1158 im[i * mz + j] = hglj(i, zp, zglj, nz, alpha, beta);
1159 }
1160 }
1161
1162 return;
1163}
Polylib::hglj
double hglj(const int i, const double z, const double *zglj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Lobatto-Jacobi points zgrj ...
Definition: Polylib.cpp:1025

References Nektar::LibUtilities::beta, and hglj().

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

◆ Imgrjm()

void Polylib::Imgrjm ( double *  im,
const double *  zgrj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Radau-Jacobi points (including z=-1) to an arbitrary distrubtion at points zm.

  • Computes the one-dimensional interpolation matrix, im, to interpolate a function from at Gauss-Radau-Jacobi distribution of nz zeros zgrj (where zgrj[0]=-1) to an arbitrary distribution of mz points zm, i.e.
    \( u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \)

Definition at line 1085 of file Polylib.cpp.

1087{
1088 double zp;
1089 int i, j;
1090
1091 for (i = 0; i < nz; i++)
1092 {
1093 for (j = 0; j < mz; j++)
1094 {
1095 zp = zm[j];
1096 im[i * mz + j] = hgrjm(i, zp, zgrj, nz, alpha, beta);
1097 }
1098 }
1099
1100 return;
1101}
Polylib::hgrjm
double hgrjm(const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at...
Definition: Polylib.cpp:951

References Nektar::LibUtilities::beta, and hgrjm().

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

◆ Imgrjp()

void Polylib::Imgrjp ( double *  im,
const double *  zgrj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Radau-Jacobi points (including z=1) to an arbitrary distrubtion at points zm.

  • Computes the one-dimensional interpolation matrix, im, to interpolate a function from at Gauss-Radau-Jacobi distribution of nz zeros zgrj (where zgrj[nz-1]=1) to an arbitrary distribution of mz points zm, i.e.
    \( u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \)

Definition at line 1116 of file Polylib.cpp.

1118{
1119 double zp;
1120 int i, j;
1121
1122 for (i = 0; i < nz; i++)
1123 {
1124 for (j = 0; j < mz; j++)
1125 {
1126 zp = zm[j];
1127 im[i * mz + j] = hgrjp(i, zp, zgrj, nz, alpha, beta);
1128 }
1129 }
1130
1131 return;
1132}
Polylib::hgrjp
double hgrjp(const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at...
Definition: Polylib.cpp:988

References Nektar::LibUtilities::beta, and hgrjp().

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

◆ jacobd()

void Polylib::jacobd ( const int  np,
const double *  z,
double *  polyd,
const int  n,
const double  alpha,
const double  beta 
)

Calculate the derivative of Jacobi polynomials.

  • Generates a vector poly of values of the derivative of the n th order Jacobi polynomial \( P^(\alpha,\beta)_n(z)\) at the np points z.
  • To do this we have used the relation
    \( \frac{d}{dz} P^{\alpha,\beta}_n(z) = \frac{1}{2} (\alpha + \beta + n + 1) P^{\alpha,\beta}_n(z) \)
  • This formulation is valid for \( -1 \leq z \leq 1 \)

Definition at line 1378 of file Polylib.cpp.

1380{
1381 int i;
1382 double one = 1.0;
1383 if (n == 0)
1384 {
1385 for (i = 0; i < np; ++i)
1386 {
1387 polyd[i] = 0.0;
1388 }
1389 }
1390 else
1391 {
1392 // jacobf(np,z,polyd,n-1,alpha+one,beta+one);
1393 jacobfd(np, z, polyd, nullptr, n - 1, alpha + one, beta + one);
1394 for (i = 0; i < np; ++i)
1395 {
1396 polyd[i] *= 0.5 * (alpha + beta + (double)n + one);
1397 }
1398 }
1399 return;
1400}
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:1248

References Nektar::LibUtilities::beta, jacobfd(), and Nektar::UnitTests::z().

Referenced by Dgj(), Dglj(), Dgrjm(), Dgrjp(), Nektar::SolverUtils::AdvectionFR::SetupCFunctions(), Nektar::SolverUtils::DiffusionLFR::SetupCFunctions(), Nektar::SolverUtils::DiffusionLFRNS::SetupCFunctions(), 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(), and zwgj().

◆ jacobfd()

void Polylib::jacobfd ( const int  np,
const double *  z,
double *  poly_in,
double *  polyd,
const int  n,
const double  alpha,
const double  beta 
)

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

  • This function returns the vectors poly_in and poly_d containing the value of the \( n^th \) order Jacobi polynomial \( P^{\alpha,\beta}_n(z) \alpha > -1, \beta > -1 \) and its derivative at the np points in z[i]
  • If poly_in = NULL then only calculate derivatice
  • If polyd = NULL then only calculate polynomial
  • To calculate the polynomial this routine uses the recursion relationship (see appendix A ref [4]) : \( \begin{array}{rcl} P^{\alpha,\beta}_0(z) &=& 1 \\ P^{\alpha,\beta}_1(z) &=& \frac{1}{2} [ \alpha-\beta+(\alpha+\beta+2)z] \\ a^1_n P^{\alpha,\beta}_{n+1}(z) &=& (a^2_n + a^3_n z) P^{\alpha,\beta}_n(z) - a^4_n P^{\alpha,\beta}_{n-1}(z) \\ a^1_n &=& 2(n+1)(n+\alpha + \beta + 1)(2n + \alpha + \beta) \\ a^2_n &=& (2n + \alpha + \beta + 1)(\alpha^2 - \beta^2) \\ a^3_n &=& (2n + \alpha + \beta)(2n + \alpha + \beta + 1) (2n + \alpha + \beta + 2) \\ a^4_n &=& 2(n+\alpha)(n+\beta)(2n + \alpha + \beta + 2) \end{array} \)
  • To calculate the derivative of the polynomial this routine uses the relationship (see appendix A ref [4]) : \( \begin{array}{rcl} b^1_n(z)\frac{d}{dz} P^{\alpha,\beta}_n(z)&=&b^2_n(z)P^{\alpha,\beta}_n(z) + b^3_n(z) P^{\alpha,\beta}_{n-1}(z) \hspace{2.2cm} \\ b^1_n(z) &=& (2n+\alpha + \beta)(1-z^2) \\ b^2_n(z) &=& n[\alpha - \beta - (2n+\alpha + \beta)z]\\ b^3_n(z) &=& 2(n+\alpha)(n+\beta) \end{array} \)
  • Note the derivative from this routine is only valid for -1 < z < 1.

Definition at line 1248 of file Polylib.cpp.

1250{
1251 int i;
1252 double zero = 0.0, one = 1.0, two = 2.0;
1253
1254 if (!np)
1255 {
1256 return;
1257 }
1258
1259 if (n == 0)
1260 {
1261 if (poly_in)
1262 {
1263 for (i = 0; i < np; ++i)
1264 {
1265 poly_in[i] = one;
1266 }
1267 }
1268 if (polyd)
1269 {
1270 for (i = 0; i < np; ++i)
1271 {
1272 polyd[i] = zero;
1273 }
1274 }
1275 }
1276 else if (n == 1)
1277 {
1278 if (poly_in)
1279 {
1280 for (i = 0; i < np; ++i)
1281 {
1282 poly_in[i] = 0.5 * (alpha - beta + (alpha + beta + two) * z[i]);
1283 }
1284 }
1285 if (polyd)
1286 {
1287 for (i = 0; i < np; ++i)
1288 {
1289 polyd[i] = 0.5 * (alpha + beta + two);
1290 }
1291 }
1292 }
1293 else
1294 {
1295 int k;
1296 double a1, a2, a3, a4;
1297 double two = 2.0, apb = alpha + beta;
1298 double *poly, *polyn1, *polyn2;
1299
1300 if (poly_in)
1301 { // switch for case of no poynomial function return
1302 polyn1 = (double *)malloc(2 * np * sizeof(double));
1303 polyn2 = polyn1 + np;
1304 poly = poly_in;
1305 }
1306 else
1307 {
1308 polyn1 = (double *)malloc(3 * np * sizeof(double));
1309 polyn2 = polyn1 + np;
1310 poly = polyn2 + np;
1311 }
1312
1313 for (i = 0; i < np; ++i)
1314 {
1315 polyn2[i] = one;
1316 polyn1[i] = 0.5 * (alpha - beta + (alpha + beta + two) * z[i]);
1317 }
1318
1319 for (k = 2; k <= n; ++k)
1320 {
1321 a1 = two * k * (k + apb) * (two * k + apb - two);
1322 a2 = (two * k + apb - one) * (alpha * alpha - beta * beta);
1323 a3 =
1324 (two * k + apb - two) * (two * k + apb - one) * (two * k + apb);
1325 a4 = two * (k + alpha - one) * (k + beta - one) * (two * k + apb);
1326
1327 a2 /= a1;
1328 a3 /= a1;
1329 a4 /= a1;
1330
1331 for (i = 0; i < np; ++i)
1332 {
1333 poly[i] = (a2 + a3 * z[i]) * polyn1[i] - a4 * polyn2[i];
1334 polyn2[i] = polyn1[i];
1335 polyn1[i] = poly[i];
1336 }
1337 }
1338
1339 if (polyd)
1340 {
1341 a1 = n * (alpha - beta);
1342 a2 = n * (two * n + alpha + beta);
1343 a3 = two * (n + alpha) * (n + beta);
1344 a4 = (two * n + alpha + beta);
1345 a1 /= a4;
1346 a2 /= a4;
1347 a3 /= a4;
1348
1349 // note polyn2 points to polyn1 at end of poly iterations
1350 for (i = 0; i < np; ++i)
1351 {
1352 polyd[i] = (a1 - a2 * z[i]) * poly[i] + a3 * polyn2[i];
1353 polyd[i] /= (one - z[i] * z[i]);
1354 }
1355 }
1356
1357 free(polyn1);
1358 }
1359
1360 return;
1361}

References Nektar::LibUtilities::beta, and Nektar::UnitTests::z().

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

◆ Jacobz()

static void Polylib::Jacobz ( const int  n,
double *  z,
const double  alpha,
const double  beta 
)
static

Calculate the n zeros, z, of the Jacobi polynomial, i.e. \( P_n^{\alpha,\beta}(z) = 0 \).

This routine is only value for \(( \alpha > -1, \beta > -1)\) and uses polynomial deflation in a Newton iteration

Definition at line 1597 of file Polylib.cpp.

1599{
1600 int i, j, k;
1601 double dth = M_PI / (2.0 * (double)n);
1602 double poly, pder, rlast = 0.0;
1603 double sum, delr, r;
1604 double one = 1.0, two = 2.0;
1605
1606 if (!n)
1607 {
1608 return;
1609 }
1610
1611 for (k = 0; k < n; ++k)
1612 {
1613 r = -cos((two * (double)k + one) * dth);
1614 if (k)
1615 {
1616 r = 0.5 * (r + rlast);
1617 }
1618
1619 for (j = 1; j < STOP; ++j)
1620 {
1621 jacobfd(1, &r, &poly, &pder, n, alpha, beta);
1622
1623 for (i = 0, sum = 0.0; i < k; ++i)
1624 {
1625 sum += one / (r - z[i]);
1626 }
1627
1628 delr = -poly / (pder - sum * poly);
1629 r += delr;
1630 if (fabs(delr) < EPS)
1631 {
1632 break;
1633 }
1634 }
1635 z[k] = r;
1636 rlast = r;
1637 }
1638 return;
1639}
STOP
#define STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
Definition: Polylib.cpp:43

References Nektar::LibUtilities::beta, EPS, jacobfd(), STOP, and Nektar::UnitTests::z().

◆ JacZeros()

void Polylib::JacZeros ( const int  n,
double *  a,
double *  b,
const double  alpha,
const double  beta 
)

Zero and Weight determination through the eigenvalues and eigenvectors of a tridiagonal matrix from the three term recurrence relationship.

Set up a symmetric tridiagonal matrix

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

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

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

where \( j=n+1\) and \(p_j(z)\) are the Jacobi (normalized) orthogonal polynomials \( \alpha,\beta > -1\)( integer values and halves). Since the polynomials are orthonormalized, the tridiagonal matrix is guaranteed to be symmetric. The eigenvalues of this matrix are the zeros of the Jacobi polynomial.

Definition at line 1665 of file Polylib.cpp.

1667{
1668
1669 int i, j;
1670 RecCoeff(n, a, b, alpha, beta);
1671
1672 double **z = new double *[n];
1673 for (i = 0; i < n; i++)
1674 {
1675 z[i] = new double[n];
1676 for (j = 0; j < n; j++)
1677 {
1678 z[i][j] = 0.0;
1679 }
1680 }
1681 for (i = 0; i < n; i++)
1682 {
1683 z[i][i] = 1.0;
1684 }
1685
1686 // find eigenvalues and eigenvectors
1687 TriQL(n, a, b, z);
1688
1689 delete[] z;
1690 return;
1691}
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and b of the orthogonal polynomials.
Definition: Polylib.cpp:1697
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:1758

References Nektar::LibUtilities::beta, RecCoeff(), TriQL(), and Nektar::UnitTests::z().

◆ JKMatrix()

void Polylib::JKMatrix ( int  n,
double *  a,
double *  b 
)

Calcualtes the Jacobi-kronrod matrix by determining the a and coefficients.

The first 3n+1 coefficients are already known

For more information refer to: "Dirk P. Laurie, Calcualtion of Gauss-Kronrod quadrature rules"

Definition at line 1875 of file Polylib.cpp.

1876{
1877 int i, j, k, m;
1878 // Working storage
1879 int size = (int)floor(n / 2.0) + 2;
1880 double *s = new double[size];
1881 double *t = new double[size];
1882
1883 // Initialize s and t to zero
1884 for (i = 0; i < size; i++)
1885 {
1886 s[i] = 0.0;
1887 t[i] = 0.0;
1888 }
1889
1890 t[1] = b[n + 1];
1891 for (m = 0; m <= n - 2; m++)
1892 {
1893 double u = 0.0;
1894 for (k = (int)floor((m + 1) / 2.0); k >= 0; k--)
1895 {
1896 int l = m - k;
1897 u = u + (a[k + n + 1] - a[l]) * t[k + 1] + b[k + n + 1] * s[k] -
1898 b[l] * s[k + 1];
1899 s[k + 1] = u;
1900 }
1901
1902 // Swap the contents of s and t
1903 double *hold = s;
1904 s = t;
1905 t = hold;
1906 }
1907
1908 for (j = (int)floor(n / 2.0); j >= 0; j--)
1909 {
1910 s[j + 1] = s[j];
1911 }
1912
1913 for (m = n - 1; m <= 2 * n - 3; m++)
1914 {
1915 double u = 0;
1916 for (k = m + 1 - n; k <= floor((m - 1) / 2.0); k++)
1917 {
1918 int l = m - k;
1919 j = n - 1 - l;
1920 u = u - (a[k + n + 1] - a[l]) * t[j + 1] - b[k + n + 1] * s[j + 1] +
1921 b[l] * s[j + 2];
1922 s[j + 1] = u;
1923 }
1924
1925 if (m % 2 == 0)
1926 {
1927 k = m / 2;
1928 a[k + n + 1] =
1929 a[k] + (s[j + 1] - b[k + n + 1] * s[j + 2]) / t[j + 2];
1930 }
1931 else
1932 {
1933 k = (m + 1) / 2;
1934 b[k + n + 1] = s[j + 1] / s[j + 2];
1935 }
1936
1937 // Swap the contents of s and t
1938 double *hold = s;
1939 s = t;
1940 t = hold;
1941 }
1942
1943 a[2 * n] = a[n - 1] - b[2 * n] * s[1] / t[1];
1944}

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

◆ laginterp()

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

Definition at line 85 of file Polylib.cpp.

86{
87 double temp = 1.0;
88 for (int i = 0; i < np; i++)
89 {
90 if (j != i)
91 {
92 temp *= optdiff(z, zj[i]) / (zj[j] - zj[i]);
93 }
94 }
95 return temp;
96}
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:55

References optdiff(), and Nektar::UnitTests::z().

Referenced by Nektar::LibUtilities::GaussPoints::CalculateInterpMatrix(), Nektar::LibUtilities::PolyEPoints::CalculateInterpMatrix(), hgj(), hglj(), hgrjm(), hgrjp(), and Nektar::LibUtilities::PolyEPoints::v_CalculateWeights().

◆ laginterpderiv()

double Polylib::laginterpderiv ( double  z,
int  k,
const double *  zj,
int  np 
)

Definition at line 98 of file Polylib.cpp.

99{
100 double y = 0.0;
101 for (int j = 0; j < np; ++j)
102 {
103 if (j != k)
104 {
105 double tmp = 1.0;
106 for (int i = 0; i < np; ++i)
107 {
108 if (i != k)
109 {
110 if (i != j)
111 {
112 tmp *= (z - zj[i]);
113 }
114 tmp /= (zj[k] - zj[i]);
115 }
116 }
117 y += tmp;
118 }
119 }
120 return y;
121}

References Nektar::UnitTests::z().

Referenced by Nektar::LibUtilities::GaussPoints::v_CalculateDerivMatrix(), and Nektar::LibUtilities::PolyEPoints::v_CalculateDerivMatrix().

◆ optdiff()

double Polylib::optdiff ( double  xl,
double  xr 
)

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

Definition at line 55 of file Polylib.cpp.

56{
57 double m_xln, m_xrn;
58 int m_expn;
59 int m_digits = static_cast<int>(fabs(floor(log10(DBL_EPSILON))) - 1);
60
61 if (fabs(xl - xr) < 1.e-4)
62 {
63
64 m_expn = static_cast<int>(floor(log10(fabs(xl - xr))));
65 m_xln = xl * powl(10.0L, -m_expn) -
66 floor(xl * powl(10.0L,
67 -m_expn)); // substract the digits overlap part
68 m_xrn =
69 xr * powl(10.0L, -m_expn) -
70 floor(xl *
71 powl(10.0L,
72 -m_expn)); // substract the common digits overlap part
73 m_xln =
74 round(m_xln * powl(10.0L, m_digits + m_expn)); // git rid of rubbish
75 m_xrn = round(m_xrn * powl(10.0L, m_digits + m_expn));
76
77 return powl(10.0L, -m_digits) * (m_xln - m_xrn);
78 }
79 else
80 {
81 return (xl - xr);
82 }
83}

Referenced by laginterp().

◆ polycoeffs()

void Polylib::polycoeffs ( double *  c,
const double *  z,
const int  i,
const int  np 
)

Compute the coefficients of Lagrange interpolation polynomials.

Definition at line 1170 of file Polylib.cpp.

1171{
1172 int j, k, m;
1173
1174 // Compute denominator
1175 double d = 1.0;
1176 for (j = 0; j < np; j++)
1177 {
1178 if (i != j)
1179 {
1180 d *= z[i] - z[j];
1181 }
1182 c[j] = 0.0;
1183 }
1184
1185 // Compute coefficient
1186 c[0] = 1.0;
1187 m = 0;
1188 for (j = 0; j < np; j++)
1189 {
1190 if (i != j)
1191 {
1192 m += 1;
1193 c[m] = c[m - 1];
1194 for (k = m - 1; k > 0; k--)
1195 {
1196 c[k] *= -1.0 * z[j];
1197 c[k] += c[k - 1];
1198 }
1199 c[0] *= -1.0 * z[j];
1200 }
1201 }
1202 for (j = 0; j < np; j++)
1203 {
1204 c[j] /= d;
1205 }
1206}
Nektar::UnitTests::d
std::vector< double > d(NPUPPER *NPUPPER)

References Nektar::UnitTests::d(), and Nektar::UnitTests::z().

Referenced by Qg().

◆ Qg()

void Polylib::Qg ( double *  Q,
const double *  z,
const int  np 
)

Compute the Integration Matrix.

Definition at line 611 of file Polylib.cpp.

612{
613
614 if (np <= 0)
615 {
616 Q[0] = 0.0;
617 }
618 else
619 {
620 int i, j, k;
621 double *pd;
622
623 pd = (double *)malloc(np * sizeof(double));
624
625 for (i = 0; i < np; i++)
626 {
627 polycoeffs(pd, z, i, np);
628 for (j = 0; j < np; j++)
629 {
630 Q[j * np + i] = 0.0;
631 for (k = 0; k < np; k++)
632 {
633 Q[j * np + i] += pd[k] / (k + 1) * pow(z[j], k + 1);
634 }
635 }
636 }
637 free(pd);
638 }
639 return;
640}
Polylib::polycoeffs
void polycoeffs(double *c, const double *z, const int i, const int np)
Compute the coefficients of Lagrange interpolation polynomials.
Definition: Polylib.cpp:1170

References polycoeffs(), and Nektar::UnitTests::z().

Referenced by Nektar::UnitTests::BOOST_AUTO_TEST_CASE(), and Nektar::LibUtilities::TimeIntegrationSchemeSDC::v_InitializeScheme().

◆ RecCoeff()

static void Polylib::RecCoeff ( const int  n,
double *  a,
double *  b,
const double  alpha,
const double  beta 
)
static

The routine finds the recurrence coefficients a and b of the orthogonal polynomials.

Definition at line 1697 of file Polylib.cpp.

1699{
1700
1701 int i;
1702 double apb, apbi, a2b2;
1703
1704 if (!n)
1705 {
1706 return;
1707 }
1708
1709 // generate normalised terms
1710 apb = alpha + beta;
1711 apbi = 2.0 + apb;
1712
1713 b[0] = pow(2.0, apb + 1.0) * gammaF(alpha + 1.0) * gammaF(beta + 1.0) /
1714 gammaF(apbi); // MuZero
1715 a[0] = (beta - alpha) / apbi;
1716 b[1] = (4.0 * (1.0 + alpha) * (1.0 + beta) / ((apbi + 1.0) * apbi * apbi));
1717
1718 a2b2 = beta * beta - alpha * alpha;
1719
1720 for (i = 1; i < n - 1; i++)
1721 {
1722 apbi = 2.0 * (i + 1) + apb;
1723 a[i] = a2b2 / ((apbi - 2.0) * apbi);
1724 b[i + 1] = (4.0 * (i + 1) * (i + 1 + alpha) * (i + 1 + beta) *
1725 (i + 1 + apb) / ((apbi * apbi - 1) * apbi * apbi));
1726 }
1727
1728 apbi = 2.0 * n + apb;
1729 a[n - 1] = a2b2 / ((apbi - 2.0) * apbi);
1730}

References Nektar::LibUtilities::beta, and gammaF().

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

◆ TriQL()

static void Polylib::TriQL ( const int  n,
double *  d,
double *  e,
double **  z 
)
static

QL algorithm for symmetric tridiagonal matrix.

This subroutine is a translation of an algol procedure, num. math. 12, 377-383(1968) by martin and wilkinson, as modified in num. math. 15, 450(1970) by dubrulle. Handbook for auto. comp., vol.ii-linear algebra, 241-248(1971). This is a modified version from numerical recipes.

This subroutine finds the eigenvalues and first components of the eigenvectors of a symmetric tridiagonal matrix by the implicit QL method.

on input:

  • n is the order of the matrix;
  • d contains the diagonal elements of the input matrix;
  • e contains the subdiagonal elements of the input matrix in its first n-2 positions.
    • z is the n by n identity matrix

on output:

  • d contains the eigenvalues in ascending order.
  • e contains the weight values - modifications of the first component of normalised eigenvectors

Definition at line 1758 of file Polylib.cpp.

1759{
1760 int m, l, iter, i, k;
1761 double s, r, p, g, f, dd, c, b;
1762
1763 double MuZero = e[0];
1764
1765 // Renumber the elements of e
1766 for (i = 0; i < n - 1; i++)
1767 {
1768 e[i] = sqrt(e[i + 1]);
1769 }
1770 e[n - 1] = 0.0;
1771
1772 for (l = 0; l < n; l++)
1773 {
1774 iter = 0;
1775 do
1776 {
1777 for (m = l; m < n - 1; m++)
1778 {
1779 dd = fabs(d[m]) + fabs(d[m + 1]);
1780 if (fabs(e[m]) + dd == dd)
1781 {
1782 break;
1783 }
1784 }
1785 if (m != l)
1786 {
1787 if (iter++ == STOP)
1788 {
1789 fprintf(stderr, "triQL: Too many iterations in TQLI");
1790 exit(1);
1791 }
1792 g = (d[l + 1] - d[l]) / (2.0 * e[l]);
1793 r = sqrt((g * g) + 1.0);
1794 g = d[m] - d[l] + e[l] / (g + sign(r, g));
1795 s = c = 1.0;
1796 p = 0.0;
1797 for (i = m - 1; i >= l; i--)
1798 {
1799 f = s * e[i];
1800 b = c * e[i];
1801 if (fabs(f) >= fabs(g))
1802 {
1803 c = g / f;
1804 r = sqrt((c * c) + 1.0);
1805 e[i + 1] = f * r;
1806 c *= (s = 1.0 / r);
1807 }
1808 else
1809 {
1810 s = f / g;
1811 r = sqrt((s * s) + 1.0);
1812 e[i + 1] = g * r;
1813 s *= (c = 1.0 / r);
1814 }
1815 g = d[i + 1] - p;
1816 r = (d[i] - g) * s + 2.0 * c * b;
1817 p = s * r;
1818 d[i + 1] = g + p;
1819 g = c * r - b;
1820
1821 // Calculate the eigenvectors
1822 for (k = 0; k < n; k++)
1823 {
1824 f = z[k][i + 1];
1825 z[k][i + 1] = s * z[k][i] + c * f;
1826 z[k][i] = c * z[k][i] - s * f;
1827 }
1828 }
1829 d[l] = d[l] - p;
1830 e[l] = g;
1831 e[m] = 0.0;
1832 }
1833 } while (m != l);
1834 }
1835
1836 // order eigenvalues
1837 // Since we only need the first component of the eigenvectors
1838 // to calcualte the weight, we only swap the first components
1839 for (i = 0; i < n - 1; ++i)
1840 {
1841 k = i;
1842 p = d[i];
1843 for (l = i + 1; l < n; ++l)
1844 {
1845 if (d[l] < p)
1846 {
1847 k = l;
1848 p = d[l];
1849 }
1850 }
1851 d[k] = d[i];
1852 d[i] = p;
1853
1854 double temp = z[0][k];
1855 z[0][k] = z[0][i];
1856 z[0][i] = temp;
1857 }
1858
1859 // Calculate the weights
1860 for (i = 0; i < n; i++)
1861 {
1862 e[i] = MuZero * z[0][i] * z[0][i];
1863 }
1864}
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:47
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350

References Nektar::UnitTests::d(), CellMLToNektar.cellml_metadata::p, sign, tinysimd::sqrt(), STOP, and Nektar::UnitTests::z().

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

◆ zwgj()

void Polylib::zwgj ( double *  z,
double *  w,
const int  np,
const double  alpha,
const double  beta 
)

Gauss-Jacobi zeros and weights.

  • Generate np Gauss Jacobi zeros, z, and weights,w, associated with the Jacobi polynomial \( P^{\alpha,\beta}_{np}(z) \),
  • Exact for polynomials of order 2np-1 or less

Definition at line 152 of file Polylib.cpp.

154{
155 int i;
156 double fac, one = 1.0, two = 2.0, apb = alpha + beta;
157
158 jacobz(np, z, alpha, beta);
159 jacobd(np, z, w, np, alpha, beta);
160
161 fac = pow(two, apb + one) * gammaFracGammaF(np + 1, alpha, np + 1, 0.0) *
162 gammaFracGammaF(np + 1, beta, np + 1, apb);
163
164 for (i = 0; i < np; ++i)
165 {
166 w[i] = fac / (w[i] * w[i] * (one - z[i] * z[i]));
167 }
168
169 return;
170}
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:128
Nektar::UnitTests::w
std::vector< double > w(NPUPPER)

References Nektar::LibUtilities::beta, gammaFracGammaF(), jacobd(), jacobz, Nektar::UnitTests::w(), and Nektar::UnitTests::z().

Referenced by Nektar::UnitTests::BOOST_AUTO_TEST_CASE(), and Nektar::LibUtilities::GaussPoints::v_CalculatePoints().

◆ zwgk()

void Polylib::zwgk ( double *  z,
double *  w,
const int  npt,
const double  alpha,
const double  beta 
)

Gauss-Kronrod-Jacobi zeros and weights.

  • Generate npt=2*np+1 Gauss-Kronrod Jacobi zeros, z, and weights,w, associated with the Jacobi polynomial \( P^{\alpha,\beta}_{np}(z) \),
  • Exact for polynomials of order 3np+1 or less

Definition at line 317 of file Polylib.cpp.

319{
320
321 int np = (npt - 1) / 2;
322
323 int i, j;
324
325 // number of kronrod points associated with the np gauss rule
326 int kpoints = 2 * np + 1;
327
328 // Define the number of required recurrence coefficents
329 int ncoeffs = (int)floor(3.0 * (np + 1) / 2);
330
331 // Define arrays for the recurrence coefficients
332 // We will use these arrays for the Kronrod results too, hence the
333 // reason for the size of the arrays
334 double *a = new double[kpoints];
335 double *b = new double[kpoints];
336
337 // Initialize a and b to zero
338 for (i = 0; i < kpoints; i++)
339 {
340 a[i] = 0.0;
341 b[i] = 0.0;
342 }
343
344 // Call the routine to calculate the recurrence coefficients
345 RecCoeff(ncoeffs, a, b, alpha, beta);
346
347 // Call the routine to caluclate the jacobi-Kronrod matrix
348 JKMatrix(np, a, b);
349
350 // Set up the identity matrix
351 double **zmatrix = new double *[kpoints];
352 for (i = 0; i < kpoints; i++)
353 {
354 zmatrix[i] = new double[kpoints];
355 for (j = 0; j < kpoints; j++)
356 {
357 zmatrix[i][j] = 0.0;
358 }
359 }
360 for (i = 0; i < kpoints; i++)
361 {
362 zmatrix[i][i] = 1.0;
363 }
364
365 // Calculte the points and weights
366 TriQL(kpoints, a, b, zmatrix);
367
368 for (i = 0; i < kpoints; i++)
369 {
370 z[i] = a[i];
371 w[i] = b[i];
372 }
373 delete[] a;
374 delete[] b;
375 for (i = 0; i < kpoints; i++)
376 {
377 delete[] zmatrix[i];
378 }
379 delete[] zmatrix;
380}
Polylib::JKMatrix
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the a and coefficients.
Definition: Polylib.cpp:1875

References Nektar::LibUtilities::beta, JKMatrix(), RecCoeff(), TriQL(), Nektar::UnitTests::w(), and Nektar::UnitTests::z().

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

◆ zwglj()

void Polylib::zwglj ( double *  z,
double *  w,
const int  np,
const double  alpha,
const double  beta 
)

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

  • Generate np Gauss-Lobatto-Jacobi points, z, and weights, w, associated with polynomial \( (1-z)(1+z) P^{\alpha+1,\beta+1}_{np-2}(z) \)
  • Exact for polynomials of order 2np-3 or less

Definition at line 266 of file Polylib.cpp.

268{
269
270 if (np == 1)
271 {
272 z[0] = 0.0;
273 w[0] = 2.0;
274 }
275 else if (np == 2)
276 {
277 z[0] = -1.0;
278 z[1] = 1.0;
279
280 w[0] = 1.0;
281 w[1] = 1.0;
282 }
283 else
284 {
285 int i;
286 double fac, one = 1.0, apb = alpha + beta, two = 2.0;
287
288 z[0] = -one;
289 z[np - 1] = one;
290 jacobz(np - 2, z + 1, alpha + one, beta + one);
291 jacobfd(np, z, w, nullptr, np - 1, alpha, beta);
292
293 fac = pow(two, apb + 1) * gammaFracGammaF(np, alpha, np, 0.0) *
294 gammaFracGammaF(np, beta, np + 1, apb);
295 fac /= (np - 1);
296
297 for (i = 0; i < np; ++i)
298 {
299 w[i] = fac / (w[i] * w[i]);
300 }
301 w[0] *= (beta + one);
302 w[np - 1] *= (alpha + one);
303 }
304
305 return;
306}

References Nektar::LibUtilities::beta, gammaFracGammaF(), jacobfd(), jacobz, Nektar::UnitTests::w(), and Nektar::UnitTests::z().

Referenced by Nektar::UnitTests::BOOST_AUTO_TEST_CASE(), and Nektar::LibUtilities::GaussPoints::v_CalculatePoints().

◆ zwgrjm()

void Polylib::zwgrjm ( double *  z,
double *  w,
const int  np,
const double  alpha,
const double  beta 
)

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

  • Generate np Gauss-Radau-Jacobi zeros, z, and weights,w, associated with the polynomial \((1+z) P^{\alpha,\beta+1}_{np-1}(z) \).
  • Exact for polynomials of order 2np-2 or less

Definition at line 182 of file Polylib.cpp.

184{
185
186 if (np == 1)
187 {
188 z[0] = 0.0;
189 w[0] = 2.0;
190 }
191 else
192 {
193 int i;
194 double fac, one = 1.0, two = 2.0, apb = alpha + beta;
195
196 z[0] = -one;
197 jacobz(np - 1, z + 1, alpha, beta + 1);
198 jacobfd(np, z, w, nullptr, np - 1, alpha, beta);
199
200 fac = pow(two, apb) * gammaFracGammaF(np, alpha, np, 0.0) *
201 gammaFracGammaF(np, beta, np + 1, apb);
202 fac /= (beta + np);
203
204 for (i = 0; i < np; ++i)
205 {
206 w[i] = fac * (1 - z[i]) / (w[i] * w[i]);
207 }
208 w[0] *= (beta + one);
209 }
210
211 return;
212}

References Nektar::LibUtilities::beta, gammaFracGammaF(), jacobfd(), jacobz, Nektar::UnitTests::w(), and Nektar::UnitTests::z().

Referenced by Nektar::UnitTests::BOOST_AUTO_TEST_CASE(), and Nektar::LibUtilities::GaussPoints::v_CalculatePoints().

◆ zwgrjp()

void Polylib::zwgrjp ( double *  z,
double *  w,
const int  np,
const double  alpha,
const double  beta 
)

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

  • Generate np Gauss-Radau-Jacobi zeros, z, and weights,w, associated with the polynomial \((1-z) P^{\alpha+1,\beta}_{np-1}(z) \).
  • Exact for polynomials of order 2np-2 or less

Definition at line 225 of file Polylib.cpp.

227{
228
229 if (np == 1)
230 {
231 z[0] = 0.0;
232 w[0] = 2.0;
233 }
234 else
235 {
236 int i;
237 double fac, one = 1.0, two = 2.0, apb = alpha + beta;
238
239 jacobz(np - 1, z, alpha + 1, beta);
240 z[np - 1] = one;
241 jacobfd(np, z, w, nullptr, np - 1, alpha, beta);
242
243 fac = pow(two, apb) * gammaFracGammaF(np, alpha, np, 0.0) *
244 gammaFracGammaF(np, beta, np + 1, apb);
245 fac /= (alpha + np);
246
247 for (i = 0; i < np; ++i)
248 {
249 w[i] = fac * (1 + z[i]) / (w[i] * w[i]);
250 }
251 w[np - 1] *= (alpha + one);
252 }
253
254 return;
255}

References Nektar::LibUtilities::beta, gammaFracGammaF(), jacobfd(), jacobz, Nektar::UnitTests::w(), and Nektar::UnitTests::z().

Referenced by Nektar::UnitTests::BOOST_AUTO_TEST_CASE(), and Nektar::LibUtilities::GaussPoints::v_CalculatePoints().

◆ zwlk()

void Polylib::zwlk ( double *  z,
double *  w,
const int  npt,
const double  alpha,
const double  beta 
)

Gauss-Lobatto-Kronrod-Jacobi zeros and weights.

  • Generate npt=2*np-1 Lobatto-Kronrod Jacobi zeros, z, and weights,w, associated with the Jacobi polynomial \( P^{\alpha,\beta}_{np}(z) \),

Definition at line 495 of file Polylib.cpp.

497{
498
499 int np = (npt + 1) / 2;
500
501 if (np < 4)
502 {
503 fprintf(stderr, "too few points in formula\n");
504 return;
505 }
506
507 double endl = -1;
508 double endr = 1;
509 int i, j;
510
511 // number of kronrod points associated with the np gauss rule
512 int kpoints = 2 * np - 1;
513
514 // Define the number of required recurrence coefficents
515 int ncoeffs = (int)ceil(3.0 * np / 2) - 1;
516
517 // Define arrays for the recurrence coefficients
518 double *a = new double[ncoeffs + 1];
519 double *b = new double[ncoeffs + 1];
520
521 // Initialize a and b to zero
522 for (i = 0; i < ncoeffs + 1; i++)
523 {
524 a[i] = 0.0;
525 b[i] = 0.0;
526 }
527
528 // Call the routine to calculate the recurrence coefficients
529 RecCoeff(ncoeffs, a, b, alpha, beta);
530
531 double *a0 = new double[ncoeffs];
532 double *b0 = new double[ncoeffs];
533
534 chri1(ncoeffs, a, b, a0, b0, endl);
535
536 double *a1 = new double[ncoeffs - 1];
537 double *b1 = new double[ncoeffs - 1];
538
539 chri1(ncoeffs - 1, a0, b0, a1, b1, endr);
540
541 double s = b1[0] / fabs(b1[0]);
542 b1[0] = s * b1[0];
543
544 // Finding the 2*np-1 gauss-kronrod points
545 double *z1 = new double[2 * np - 3];
546 double *w1 = new double[2 * np - 3];
547 for (i = 0; i < ncoeffs; i++)
548 {
549 z1[i] = a1[i];
550 w1[i] = b1[i];
551 }
552 JKMatrix(np - 2, z1, w1);
553 // Set up the identity matrix
554 double **zmatrix = new double *[2 * np - 3];
555 for (i = 0; i < 2 * np - 3; i++)
556 {
557 zmatrix[i] = new double[2 * np - 3];
558 for (j = 0; j < 2 * np - 3; j++)
559 {
560 zmatrix[i][j] = 0.0;
561 }
562 }
563 for (i = 0; i < 2 * np - 3; i++)
564 {
565 zmatrix[i][i] = 1.0;
566 }
567
568 // Calculate the points and weights
569 TriQL(2 * np - 3, z1, w1, zmatrix);
570
571 double sumW1 = 0.0;
572 double sumW1Z1 = 0.0;
573 for (i = 0; i < 2 * np - 3; i++)
574 {
575 w1[i] = s * w1[i] / (z1[i] - endl) / (z1[i] - endr);
576 sumW1 += w1[i];
577 sumW1Z1 += z1[i] * w1[i];
578 }
579
580 double c0 = b[0] - sumW1;
581 double c1 = a[0] * b[0] - sumW1Z1;
582
583 z[0] = endl;
584 z[2 * np - 2] = endr;
585 w[0] = (c0 * endr - c1) / (endr - endl);
586 w[2 * np - 2] = (c1 - c0 * endl) / (endr - endl);
587
588 for (i = 1; i < kpoints - 1; i++)
589 {
590 z[i] = z1[i - 1];
591 w[i] = w1[i - 1];
592 }
593 delete[] a;
594 delete[] b;
595 delete[] a0;
596 delete[] b0;
597 delete[] a1;
598 delete[] b1;
599 delete[] z1;
600 delete[] w1;
601 for (i = 0; i < 2 * np - 3; i++)
602 {
603 delete[] zmatrix[i];
604 }
605 delete[] zmatrix;
606}
Polylib::chri1
void chri1(int, double *, double *, double *, double *, double)
Given a weight function through the first n+1 coefficients a and b of its orthogonal polynomials thi...
Definition: Polylib.cpp:1956

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

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

◆ zwrk()

void Polylib::zwrk ( double *  z,
double *  w,
const int  npt,
const double  alpha,
const double  beta 
)

Gauss-Radau-Kronrod-Jacobi zeros and weights.

  • Generate npt=2*np Radau-Kronrod Jacobi zeros, z, and weights,w, associated with the Jacobi polynomial \( P^{\alpha,\beta}_{np}(z) \),

Definition at line 389 of file Polylib.cpp.

391{
392
393 int np = npt / 2;
394
395 if (np < 2)
396 {
397 fprintf(stderr, "too few points in formula\n");
398 return;
399 }
400
401 double end0 = -1;
402
403 int i, j;
404
405 // number of kronrod points associated with the np gauss rule
406 int kpoints = 2 * np;
407
408 // Define the number of required recurrence coefficents
409 int ncoeffs = (int)ceil(3.0 * np / 2);
410
411 // Define arrays for the recurrence coefficients
412 double *a = new double[ncoeffs + 1];
413 double *b = new double[ncoeffs + 1];
414
415 // Initialize a and b to zero
416 for (i = 0; i < ncoeffs + 1; i++)
417 {
418 a[i] = 0.0;
419 b[i] = 0.0;
420 }
421
422 // Call the routine to calculate the recurrence coefficients
423 RecCoeff(ncoeffs, a, b, alpha, beta);
424
425 double *a0 = new double[ncoeffs];
426 double *b0 = new double[ncoeffs];
427
428 chri1(ncoeffs, a, b, a0, b0, end0);
429
430 double s = b0[0] / fabs(b0[0]);
431 b0[0] = s * b0[0];
432
433 // Finding the 2*np-1 gauss-kronrod points
434 double *z1 = new double[2 * np - 1];
435 double *w1 = new double[2 * np - 1];
436 for (i = 0; i < ncoeffs; i++)
437 {
438 z1[i] = a0[i];
439 w1[i] = b0[i];
440 }
441 JKMatrix(np - 1, z1, w1);
442 // Set up the identity matrix
443 double **zmatrix = new double *[2 * np - 1];
444 for (i = 0; i < 2 * np - 1; i++)
445 {
446 zmatrix[i] = new double[2 * np - 1];
447 for (j = 0; j < 2 * np - 1; j++)
448 {
449 zmatrix[i][j] = 0.0;
450 }
451 }
452 for (i = 0; i < 2 * np - 1; i++)
453 {
454 zmatrix[i][i] = 1.0;
455 }
456
457 // Calculate the points and weights
458 TriQL(2 * np - 1, z1, w1, zmatrix);
459
460 double sumW1 = 0.0;
461 for (i = 0; i < 2 * np - 1; i++)
462 {
463 w1[i] = s * w1[i] / (z1[i] - end0);
464 sumW1 += w1[i];
465 }
466
467 z[0] = end0;
468 w[0] = b[0] - sumW1;
469 for (i = 1; i < kpoints; i++)
470 {
471 z[i] = z1[i - 1];
472 w[i] = w1[i - 1];
473 }
474
475 delete[] a;
476 delete[] b;
477 delete[] a0;
478 delete[] b0;
479 delete[] z1;
480 delete[] w1;
481 for (i = 0; i < 2 * np - 1; i++)
482 {
483 delete[] zmatrix[i];
484 }
485 delete[] zmatrix;
486}

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

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

Generated on Wed Jan 31 2024 09:36:07 for Nektar++ by doxygen 1.9.4