PolynomialQuadrature Namespace Reference

Polynomials and quadratures based on defined distributions for Polynomial Chaos. More...

Classes
class	ClenshawCurtisGrid
	Clenshaw-Curtis sparse grid. More...
class	Hermite
	Normal distributions use Hermite polynomials. More...
class	Legendre
	Uniform distributions use Legendre polynomials. More...
class	Polynomial
	General polynomial class with function for evaluating a polynomial of a given order at a given point. More...
class	Quadrature
	General multidimensional quadrature class. More...
class	SmolyakGrid
	Smolyak sparse grid. More...
class	TensorGrid
	Full tensor product of 1D quadratures. More...

Functions
std::unique_ptr< const Polynomial >	makePolynomial (const Distribution *dist)
Real	legendre (const unsigned int order, const Real x, const Real lower_bound=-1.0, const Real upper_bound=1.0)
	Legendre polynomial of specified order.
Real	hermite (const unsigned int order, const Real x, const Real mu=0.0, const Real sig=1.0)
	Hermite polynomial of specified order.
void	gauss_legendre (const unsigned int order, std::vector< Real > &points, std::vector< Real > &weights, const Real lower_bound, const Real upper_bound)
	Generalized formula for any polynomial order.
void	gauss_hermite (const unsigned int order, std::vector< Real > &points, std::vector< Real > &weights, const Real mu, const Real sig)
	Generalized formula for any polynomial order.
void	clenshaw_curtis (const unsigned int order, std::vector< Real > &points, std::vector< Real > &weights)

Detailed Description

Polynomials and quadratures based on defined distributions for Polynomial Chaos.

Function Documentation

clenshaw_curtis()

void PolynomialQuadrature::clenshaw_curtis	(	const unsigned int	order,
		std::vector< Real > &	points,
		std::vector< Real > &	weights
	)

Definition at line 395 of file PolynomialQuadrature.C.

{
  // Number of points needed
  unsigned int N = order + (order % 2);
  points.resize(N + 1);
  weights.resize(N + 1);
 
  if (N == 0)
  {
    points[0] = 0;
    weights[0] = 1;
    return;
  }
 
  std::vector<Real> dk(N / 2 + 1);
  for (unsigned int k = 0; k <= (N / 2); ++k)
    dk[k] = ((k == 0 || k == (N / 2)) ? 1.0 : 2.0) / (1.0 - 4.0 * (Real)k * (Real)k);
 
  for (unsigned int n = 0; n <= (N / 2); ++n)
  {
    Real theta = (Real)n * M_PI / ((Real)N);
    points[n] = -std::cos(theta);
    for (unsigned int k = 0; k <= (N / 2); ++k)
    {
      Real Dnk =
          ((n == 0 || n == (N / 2)) ? 0.5 : 1.0) * std::cos((Real)k * theta * 2.0) / ((Real)N);
      weights[n] += Dnk * dk[k];
    }
  }
 
  for (unsigned int n = 0; n < (N / 2); ++n)
  {
    points[N - n] = -points[n];
    weights[N - n] = weights[n];
  }
  weights[N / 2] *= 2.0;
}

Referenced by PolynomialQuadrature::Legendre::clenshawQuadrature(), and TEST().

gauss_hermite()

void PolynomialQuadrature::gauss_hermite	(	const unsigned int	order,
		std::vector< Real > &	points,
		std::vector< Real > &	weights,
		const Real	mu,
		const Real	sig
	)

Generalized formula for any polynomial order.

Resulting number of points is then: N = order + 1 The sum of weights is sqrt(2*\pi) Uses the Golub-Welsch algorithm:

• Golub, Gene & Welsch, John `‘Calculation of Gauss Quadrature Rules’'. https://web.stanford.edu/class/cme335/S0025-5718-69-99647-1.pdf

Definition at line 365 of file PolynomialQuadrature.C.

{
  // Number of points needed
  unsigned int n = order + 1;
  points.resize(n);
  weights.resize(n);
 
  DenseMatrix<Real> mat(n, n);
  DenseVector<Real> lambda(n);
  DenseVector<Real> lambdai(n);
  DenseMatrix<Real> vec(n, n);
  for (unsigned int i = 1; i < n; ++i)
  {
    mat(i, i - 1) = std::sqrt(static_cast<Real>(i));
    mat(i - 1, i) = mat(i, i - 1);
  }
  mat.evd_right(lambda, lambdai, vec);
 
  for (unsigned int i = 0; i < n; ++i)
  {
    points[i] = mu + lambda(i) * sig;
    weights[i] = vec(0, i) * vec(0, i);
  }
}

Referenced by PolynomialQuadrature::Hermite::gaussQuadrature(), and TEST().

gauss_legendre()

void PolynomialQuadrature::gauss_legendre	(	const unsigned int	order,
		std::vector< Real > &	points,
		std::vector< Real > &	weights,
		const Real	lower_bound,
		const Real	upper_bound
	)

Generalized formula for any polynomial order.

Resulting number of points is then: N = order + 1 The sum of weights is 2 Uses the Golub-Welsch algorithm:

• Golub, Gene & Welsch, John `‘Calculation of Gauss Quadrature Rules’'. https://web.stanford.edu/class/cme335/S0025-5718-69-99647-1.pdf

Definition at line 333 of file PolynomialQuadrature.C.

{
  unsigned int n = order + 1;
  points.resize(n);
  weights.resize(n);
 
  DenseMatrix<Real> mat(n, n);
  DenseVector<Real> lambda(n);
  DenseVector<Real> lambdai(n);
  DenseMatrix<Real> vec(n, n);
  for (unsigned int i = 1; i < n; ++i)
  {
    Real ri = i;
    mat(i, i - 1) = ri / std::sqrt(((2. * ri - 1.) * (2. * ri + 1.)));
    mat(i - 1, i) = mat(i, i - 1);
  }
  mat.evd_right(lambda, lambdai, vec);
 
  Real dx = (upper_bound - lower_bound) / 2.0;
  Real xav = (upper_bound + lower_bound) / 2.0;
  for (unsigned int i = 0; i < n; ++i)
  {
    points[i] = lambda(i) * dx + xav;
    weights[i] = vec(0, i) * vec(0, i);
  }
}

Referenced by PolynomialQuadrature::Legendre::gaussQuadrature(), and TEST().

hermite()

Real PolynomialQuadrature::hermite	(	const unsigned int	order,
		const Real	x,
		const Real	mu = 0.0,
		const Real	sig = 1.0
	)

Hermite polynomial of specified order.

Uses boost if available.

Definition at line 297 of file PolynomialQuadrature.C.

{
  Real xref = (x - mu) / sig;
#ifdef LIBMESH_HAVE_EXTERNAL_BOOST
  // Using Hermite polynomials from boost library (if available)
  // https://www.boost.org/doc/libs/1_46_1/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_poly/hermite.html
  // Need to do some modification since boost does physicists hermite polynomials:
  // H_n^prob(x) = 2^(-n/2)H_n^phys(x / sqrt(2))
  xref /= M_SQRT2; // 1 / sqrt(2)
  Real val = boost::math::hermite(order, xref);
  val /= pow(M_SQRT2, order); // 2^(-order / 2)
  return val;
#else
  // Using explicit expression of polynomial coefficients:
  // H_n(x) = n!\sum_{m=1}^{floor(n/2)} (-1)^m / (m!(n-2m)!) * x^(n-2m) / 2^m
  // https://en.wikipedia.org/wiki/Hermite_polynomials
  if (order < 13)
  {
    Real val = 0;
    for (unsigned int m = 0; m <= (order % 2 == 0 ? order / 2 : (order - 1) / 2); ++m)
    {
      Real sgn = (m % 2 == 0 ? 1.0 : -1.0);
      Real coeff =
          1.0 / Real(Utility::factorial(m) * Utility::factorial(order - 2 * m)) / pow(2.0, m);
      unsigned int ord = order - 2 * m;
      val += sgn * coeff * pow(xref, ord);
    }
    return val * Utility::factorial(order);
  }
  else
    return xref * hermite(order - 1, xref) -
           (static_cast<Real>(order) - 1.0) * hermite(order - 2, xref);
#endif
}

Referenced by PolynomialQuadrature::Hermite::compute(), PolynomialQuadrature::Hermite::computeDerivative(), hermite(), and TEST().

legendre()

Real PolynomialQuadrature::legendre	(	const unsigned int	order,
		const Real	x,
		const Real	lower_bound = -1.0,
		const Real	upper_bound = 1.0
	)

Legendre polynomial of specified order.

Uses boost if available

Definition at line 206 of file PolynomialQuadrature.C.

{
  Real xref = 2 / (upper_bound - lower_bound) * (x - (upper_bound + lower_bound) / 2);
#ifdef LIBMESH_HAVE_EXTERNAL_BOOST
  // Using Legendre polynomials from boost library (if available)
  // https://www.boost.org/doc/libs/1_46_1/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_poly/legendre.html
  return boost::math::legendre_p(order, xref);
#else
  // Using explicit expression of polynomial coefficients:
  // P_n(x) = 1/2^n\sum_{k=1}^{floor(n/2)} (-1)^k * nchoosek(n, k) * nchoosek(2n-2k, n) * x^(n-2k)
  // https://en.wikipedia.org/wiki/Legendre_polynomials
  if (order < 16)
  {
    Real val = 0;
    for (unsigned int k = 0; k <= (order % 2 == 0 ? order / 2 : (order - 1) / 2); ++k)
    {
      Real coeff =
          Real(Utility::binomial(order, k)) * Real(Utility::binomial(2 * order - 2 * k, order));
      Real sgn = (k % 2 == 0 ? 1.0 : -1.0);
      unsigned int ord = order - 2 * k;
      val += sgn * coeff * pow(xref, ord);
    }
    return val / pow(2.0, order);
  }
  else
  {
    Real ord = order;
    return ((2.0 * ord - 1.0) * xref * legendre(order - 1, xref) -
            (ord - 1.0) * legendre(order - 2, xref)) /
           ord;
  }
#endif
}

Referenced by PolynomialQuadrature::Legendre::compute(), PolynomialQuadrature::Legendre::computeDerivativeRef(), legendre(), and TEST().

makePolynomial()

std::unique_ptr< const Polynomial > PolynomialQuadrature::makePolynomial

(

const Distribution *

dist

)

Definition at line 34 of file PolynomialQuadrature.C.

{
  const Uniform * u_dist = dynamic_cast<const Uniform *>(dist);
  if (u_dist)
    return std::make_unique<const Legendre>(dist->getParam<Real>("lower_bound"),
                                            dist->getParam<Real>("upper_bound"));
 
  const Normal * n_dist = dynamic_cast<const Normal *>(dist);
  if (n_dist)
    return std::make_unique<const Hermite>(dist->getParam<Real>("mean"),
                                           dist->getParam<Real>("standard_deviation"));
 
  ::mooseError("Polynomials for '", dist->type(), "' distributions have not been implemented.");
  return nullptr;
}

Referenced by PolynomialChaosTrainer::PolynomialChaosTrainer(), and QuadratureSampler::QuadratureSampler().