From 6af60d8a90bb595aed7f601e814019361f49a774 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 6 Apr 2018 11:07:26 +0200 Subject: [PATCH] Use Jacobi roots to simplify quadrature implementation. --- source/base/quadrature_lib.cc | 236 +++------------------------------- 1 file changed, 21 insertions(+), 215 deletions(-) diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index bf7d86d59c..b13d48904a 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -15,6 +15,7 @@ #include #include +#include #include #include @@ -67,97 +68,20 @@ QGauss<1>::QGauss (const unsigned int n) if (n == 0) return; - const unsigned int m = (n+1)/2; - - // tolerance for the Newton - // iteration below. we need to make - // it adaptive since on some - // machines (for example PowerPC) - // long double is the same as - // double -- in that case we can - // only get to a certain multiple - // of the accuracy of double there, - // while on other machines we'd - // like to go further down - // - // the situation is complicated by - // the fact that even if long - // double exists and is described - // by std::numeric_limits, we may - // not actually get the additional - // precision. One case where this - // happens is on x86, where one can - // set hardware flags that disable - // long double precision even for - // long double variables. these - // flags are not usually set, but - // for example matlab sets them and - // this then breaks deal.II code - // that is run as a subroutine to - // matlab... - // - // a similar situation exists, btw, - // when running programs under - // valgrind up to and including at - // least version 3.3: valgrind's - // emulator only supports 64 bit - // arithmetic, even for 80 bit long - // doubles. - const long double - long_double_eps = static_cast(std::numeric_limits::epsilon()), - double_eps = static_cast(std::numeric_limits::epsilon()); - - // now check whether long double is more - // accurate than double, and set - // tolerances accordingly. generate a one - // that really is generated at run-time - // and is not optimized away by the - // compiler. that makes sure that the - // tolerance is set at run-time with the - // current behavior, not at compile-time - // (not doing so leads to trouble with - // valgrind for example). - volatile long double runtime_one = 1.0; - const long double tolerance - = (runtime_one + long_double_eps != runtime_one - ? - std::max (double_eps / 100, long_double_eps * 5) - : - double_eps * 5 - ); - - - for (unsigned int i=1; i<=m; ++i) - { - long double z = std::cos(numbers::PI * (i-.25)/(n+.5)); - - long double pp; - long double p1; - - // Newton iteration - do - { - // compute L_n (z) - p1 = 1.; - long double p2 = 0.; - for (unsigned int j=0; j tolerance); - - double x = .5*z; - this->quadrature_points[i-1] = Point<1>(.5-x); - this->quadrature_points[n-i] = Point<1>(.5+x); + std::vector points + = Polynomials::jacobi_polynomial_roots(n, 0, 0); - double w = 1./((1.-z*z)*pp*pp); - this->weights[i-1] = w; - this->weights[n-i] = w; + for (unsigned int i=0; i<(points.size()+1)/2; ++i) + { + this->quadrature_points[i][0] = points[i]; + this->quadrature_points[n-i-1][0] = 1.-points[i]; + + // derivative of Jacobi polynomial + const long double pp = 0.5*(n + 1)*Polynomials::jacobi_polynomial_value(n-1, 1, 1, points[i]); + const long double x = -1. + 2.*points[i]; + const double w = 1./((1.-x*x)*pp*pp); + this->weights[i] = w; + this->weights[n-i-1] = w; } } @@ -165,42 +89,6 @@ namespace internal { namespace QGaussLobatto { - /* - * Evaluate a Jacobi polynomial $ P^{\alpha, \beta}_n(x) $ specified by the - * parameters @p alpha, @p beta, @p n. Note: The Jacobi polynomials are not - * orthonormal and defined on the interval $[-1, +1]$. @p x is the point of - * evaluation. - */ - long double JacobiP(const long double x, - const int alpha, - const int beta, - const unsigned int n) - { - // the Jacobi polynomial is evaluated - // using a recursion formula. - std::vector p(n+1); - - // initial values P_0(x), P_1(x): - p[0] = 1.0L; - if (n==0) return p[0]; - p[1] = ((alpha+beta+2)*x + (alpha-beta))/2; - if (n==1) return p[1]; - - for (unsigned int i=1; i<=(n-1); ++i) - { - const int v = 2*i + alpha + beta; - const int a1 = 2*(i+1)*(i + alpha + beta + 1)*v; - const int a2 = (v + 1)*(alpha*alpha - beta*beta); - const int a3 = v*(v + 1)*(v + 2); - const int a4 = 2*(i+alpha)*(i+beta)*(v + 2); - - p[i+1] = static_cast( (a2 + a3*x)*p[i] - a4*p[i-1])/a1; - } // for - return p[n]; - } - - - /** * Evaluate the Gamma function $ \Gamma(n) = (n-1)! $. * @param n point of evaluation (integer). @@ -215,89 +103,6 @@ namespace internal - /** - * Compute Legendre-Gauss-Lobatto quadrature points in the interval $[-1, - * +1]$. They are equal to the roots of the corresponding Jacobi polynomial - * (specified by @p alpha, @p beta). @p q is the number of points. - * - * @return Vector containing nodes. - */ - std::vector - compute_quadrature_points(const unsigned int q, - const int alpha, - const int beta) - { - const unsigned int m = q-2; // no. of inner points - std::vector x(m); - - // compute quadrature points with - // a Newton algorithm. - - // Set tolerance. See class QGauss - // for detailed explanation. - const long double - long_double_eps = static_cast(std::numeric_limits::epsilon()), - double_eps = static_cast(std::numeric_limits::epsilon()); - - // check whether long double is - // more accurate than double, and - // set tolerances accordingly - volatile long double runtime_one = 1.0; - const long double tolerance - = (runtime_one + long_double_eps != runtime_one - ? - std::max (double_eps / 100, long_double_eps * 5) - : - double_eps * 5 - ); - - // The following implementation - // follows closely the one given in - // the appendix of the book by - // Karniadakis and Sherwin: - // Spectral/hp element methods for - // computational fluid dynamics - // (Oxford University Press, 2005) - - // we take the zeros of the Chebyshev - // polynomial (alpha=beta=-0.5) as - // initial values: - for (unsigned int i=0; i0) - r = (r + x[k-1])/2; - - do - { - s = 0.; - for (unsigned int i=0; i= tolerance); - - x[k] = r; - } // for - - // add boundary points: - x.insert(x.begin(), -1.L); - x.push_back(+1.L); - - return x; - } - - - /** * Compute Legendre-Gauss-Lobatto quadrature weights. The quadrature points * and weights are related to Jacobi polynomial specified by @p alpha, @p @@ -319,7 +124,7 @@ namespace internal ((q-1)*gamma(q)*gamma(alpha+beta+q+1)); for (unsigned int i=0; i::QGaussLobatto (const unsigned int n) Assert (n >= 2, ExcNotImplemented()); std::vector points - = internal::QGaussLobatto::compute_quadrature_points(n, 1, 1); + = Polynomials::jacobi_polynomial_roots(n-2, 1, 1); + points.insert(points.begin(), 0); + points.push_back(1.); std::vector w = internal::QGaussLobatto::compute_quadrature_weights(points, 0, 0); - // scale points to the interval - // [0.0, 1.0]: + // scale weights to the interval [0.0, 1.0]: for (unsigned int i=0; iquadrature_points[i] = Point<1>(0.5 + 0.5*static_cast(points[i])); - this->weights[i] = 0.5*w[i]; + this->quadrature_points[i][0] = points[i]; + this->weights[i] = 0.5*w[i]; } } -- 2.39.5