From: Martin Kronbichler Date: Tue, 3 Jan 2017 17:34:06 +0000 (+0100) Subject: Switch Legendre polynomials to stable evaluation via root representation X-Git-Tag: v8.5.0-rc1~293^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=798c52d71f269064eee9a40b1d65306315367e7f;p=dealii.git Switch Legendre polynomials to stable evaluation via root representation --- diff --git a/doc/news/changes/minor/20170103MartinKronbichler b/doc/news/changes/minor/20170103MartinKronbichler new file mode 100644 index 0000000000..d5453efcf5 --- /dev/null +++ b/doc/news/changes/minor/20170103MartinKronbichler @@ -0,0 +1,6 @@ +Improved: The representation of the two polynomial classes +Polynomials::Legendre and Polynomials::HermiteInterpolation has been changed +to the root form, which ensures stable evaluation at high degrees as opposed +to the coefficient form previously used. +
+(Martin Kronbichler, 2017/01/03) diff --git a/include/deal.II/base/polynomial.h b/include/deal.II/base/polynomial.h index 8e4f22a162..050b0cf3e7 100644 --- a/include/deal.II/base/polynomial.h +++ b/include/deal.II/base/polynomial.h @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2000 - 2016 by the deal.II authors +// Copyright (C) 2000 - 2017 by the deal.II authors // // This file is part of the deal.II library. // @@ -43,11 +43,20 @@ namespace Polynomials /** * Base class for all 1D polynomials. A polynomial is represented in this * class by its coefficients, which are set through the constructor or by - * derived classes. Evaluation of a polynomial happens through the Horner - * scheme which provides both numerical stability and a minimal number of - * numerical operations. + * derived classes. * - * @author Ralf Hartmann, Guido Kanschat, 2000, 2006 + * There are two paths for evaluation of polynomials. One is based on the + * coefficients which are evaluated through the Horner scheme which is a + * robust general-purpose scheme. An alternative and more stable evaluation + * of high-degree polynomials with roots in the unit interval is provided by + * a product in terms of the roots. This form is available for special + * polynomials such as Lagrange polynomials or Legendre polynomials and used + * with the respective constructor. To obtain this more stable evaluation + * form, the constructor with the roots in form of a Lagrange polynomial + * must be used. In case a manipulation is done that changes the roots, the + * representation is switched to the coefficient form. + * + * @author Ralf Hartmann, Guido Kanschat, 2000, 2006, Martin Kronbichler, 2011, 2017 */ template class Polynomial : public Subscriptor @@ -69,7 +78,7 @@ namespace Polynomials Polynomial (const unsigned int n); /** - * Constructor for Lagrange polynomial and its point of evaluation. The + * Constructor for a Lagrange polynomial and its point of evaluation. The * idea is to construct $\prod_{i\neq j} \frac{x-x_i}{x_j-x_i}$, where j * is the evaluation point specified as argument and the support points * contain all points (including x_j, which will internally not be @@ -87,7 +96,8 @@ namespace Polynomials * Return the value of this polynomial at the given point. * * This function uses the Horner scheme for numerical stability of the - * evaluation. + * evaluation for polynomials in the coefficient form or the product of + * terms involving the roots if that representation is used. */ number value (const number x) const; @@ -98,7 +108,8 @@ namespace Polynomials * thus determined by the size of the array passed. * * This function uses the Horner scheme for numerical stability of the - * evaluation. + * evaluation for polynomials in the coefficient form or the product of + * terms involving the roots if that representation is used. */ void value (const number x, std::vector &values) const; @@ -111,7 +122,8 @@ namespace Polynomials * space for @p n_derivatives + 1 values. * * This function uses the Horner scheme for numerical stability of the - * evaluation. + * evaluation for polynomials in the coefficient form or the product of + * terms involving the roots if that representation is used. */ void value (const number x, const unsigned int n_derivatives, @@ -361,8 +373,9 @@ namespace Polynomials /** * Legendre polynomials of arbitrary degree. Constructing a Legendre - * polynomial of degree p, the coefficients will be computed by the - * three-term recursion formula. + * polynomial of degree p, the roots will be computed by the Gauss + * formula of the respective number of points and a representation of the + * polynomial by its roots. * * @note The polynomials defined by this class differ in two aspects by what * is usually referred to as Legendre polynomials: (i) This classes defines @@ -390,38 +403,6 @@ namespace Polynomials static std::vector > generate_complete_basis (const unsigned int degree); - - private: - /** - * Coefficients for the interval $[0,1]$. - */ - static std::vector > > shifted_coefficients; - - /** - * Vector with already computed coefficients. For each degree of the - * polynomial, we keep one pointer to the list of coefficients; we do so - * rather than keeping a vector of vectors in order to simplify - * programming multithread-safe. In order to avoid memory leak, we use a - * shared_ptr in order to correctly free the memory of the vectors when - * the global destructor is called. - */ - static std::vector > > recursive_coefficients; - - /** - * Compute coefficients recursively. The coefficients are stored in a - * static data vector to be available when needed next time. Since the - * recursion is performed for the interval $[-1,1]$, the polynomials are - * shifted to $[0,1]$ by the scale and shift functions - * of Polynomial, afterwards. - */ - static void compute_coefficients (const unsigned int p); - - /** - * Get coefficients for constructor. This way, it can use the non- - * standard constructor of Polynomial. - */ - static const std::vector & - get_coefficients (const unsigned int k); }; /** @@ -576,7 +557,7 @@ namespace Polynomials * Legendre polynomials of increasing order. The implementation is * @f{align*}{ * p_0(x) &= 2x^3-3x^2+1 \\ - * p_1(x) &= -2x^2+3x^2 \\ + * p_1(x) &= -2x^3+3x^2 \\ * p_2(x) &= x^3-2x^2+x \\ * p_3(x) &= x^3-x^2 \\ * p_4(x) &= 16x^2(x-1)^2 \\ diff --git a/source/base/polynomial.cc b/source/base/polynomial.cc index e921e0dff2..546fea9d0c 100644 --- a/source/base/polynomial.cc +++ b/source/base/polynomial.cc @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2000 - 2014 by the deal.II authors +// Copyright (C) 2000 - 2017 by the deal.II authors // // This file is part of the deal.II library. // @@ -17,6 +17,7 @@ #include #include #include +#include #include #include @@ -840,177 +841,30 @@ namespace Polynomials // ------------------ class Legendre --------------- // -// Reserve space for polynomials up to degree 19. Should be sufficient -// for the start. - std::vector > > - Legendre::recursive_coefficients(20); - std::vector > > - Legendre::shifted_coefficients(20); - Legendre::Legendre (const unsigned int k) : - Polynomial (get_coefficients(k)) - {} - - - - void - Legendre::compute_coefficients (const unsigned int k_) + Polynomial (0) { - // make sure we call the - // Polynomial::shift function - // only with an argument with - // which it will not crash the - // compiler -#ifdef DEAL_II_LONG_DOUBLE_LOOP_BUG - typedef double SHIFT_TYPE; -#else - typedef long double SHIFT_TYPE; -#endif - - unsigned int k = k_; - - // first make sure that no other - // thread intercepts the - // operation of this function; - // for this, acquire the lock - // until we quit this function - Threads::Mutex::ScopedLock lock(coefficients_lock); + this->coefficients.clear(); + this->in_lagrange_product_form = true; + this->lagrange_support_points.resize(k); - // The first 2 coefficients are hard-coded - if (k==0) - k=1; - // check: does the information - // already exist? - if ((recursive_coefficients.size() < k+1) || - ((recursive_coefficients.size() >= k+1) && - (recursive_coefficients[k] == - std_cxx11::shared_ptr >()))) - // no, then generate the - // respective coefficients + // the roots of a Legendre polynomial are exactly the points in the + // Gauss-Legendre quadrature formula + if (k > 0) { - // make sure that there is enough - // space in the array for the - // coefficients, so we have to resize - // it to size k+1 - - // but it's more complicated than - // that: we call this function - // recursively, so if we simply - // resize it to k+1 here, then - // compute the coefficients for - // degree k-1 by calling this - // function recursively, then it will - // reset the size to k -- not enough - // for what we want to do below. the - // solution therefore is to only - // resize the size if we are going to - // *increase* it - if (recursive_coefficients.size() < k+1) - recursive_coefficients.resize (k+1); - - if (k<=1) - { - // create coefficients - // vectors for k=0 and k=1 - // - // allocate the respective - // amount of memory and - // later assign it to the - // coefficients array to - // make it const - std::vector *c0 = new std::vector(1); - (*c0)[0] = 1.; - - std::vector *c1 = new std::vector(2); - (*c1)[0] = 0.; - (*c1)[1] = 1.; - - // now make these arrays - // const. use shared_ptr for - // recursive_coefficients because - // that avoids a memory leak that - // would appear if we used plain - // pointers. - recursive_coefficients[0] = - std_cxx11::shared_ptr >(c0); - recursive_coefficients[1] = - std_cxx11::shared_ptr >(c1); - - // Compute polynomials - // orthogonal on [0,1] - c0 = new std::vector(*c0); - c1 = new std::vector(*c1); - - Polynomial::shift (*c0, -1.); - Polynomial::scale(*c0, 2.); - Polynomial::shift (*c1, -1.); - Polynomial::scale(*c1, 2.); - Polynomial::multiply(*c1, std::sqrt(3.)); - shifted_coefficients[0]=std_cxx11::shared_ptr >(c0); - shifted_coefficients[1]=std_cxx11::shared_ptr >(c1); - } - else - { - // for larger numbers, - // compute the coefficients - // recursively. to do so, - // we have to release the - // lock temporarily to - // allow the called - // function to acquire it - // itself - coefficients_lock.release (); - compute_coefficients(k-1); - coefficients_lock.acquire (); - - std::vector *ck = new std::vector(k+1); - - const double a = 1./(k); - const double b = a*(2*k-1); - const double c = a*(k-1); - - (*ck)[k] = b*(*recursive_coefficients[k-1])[k-1]; - (*ck)[k-1] = b*(*recursive_coefficients[k-1])[k-2]; - for (unsigned int i=1 ; i<= k-2 ; ++i) - (*ck)[i] = b*(*recursive_coefficients[k-1])[i-1] - -c*(*recursive_coefficients[k-2])[i]; - - (*ck)[0] = -c*(*recursive_coefficients[k-2])[0]; - - // finally assign the newly - // created vector to the - // const pointer in the - // coefficients array - recursive_coefficients[k] = - std_cxx11::shared_ptr >(ck); - // and compute the - // coefficients for [0,1] - ck = new std::vector(*ck); - Polynomial::shift (*ck, -1.); - Polynomial::scale(*ck, 2.); - Polynomial::multiply(*ck, std::sqrt(2.*k+1.)); - shifted_coefficients[k] = - std_cxx11::shared_ptr >(ck); - }; - }; - } - - - - const std::vector & - Legendre::get_coefficients (const unsigned int k) - { - // first make sure the coefficients - // get computed if so necessary - compute_coefficients (k); + QGauss<1> gauss(k); + for (unsigned int i=0; ilagrange_support_points[i] = gauss.get_points()[i][0]; + } - // then get a pointer to the array - // of coefficients. do that in a MT - // safe way - Threads::Mutex::ScopedLock lock (coefficients_lock); - return *shifted_coefficients[k]; + // compute the abscissa in zero of the product of monomials. The exact + // value should be sqrt(2*k+1), so set the weight to that value. + double prod = 1.; + for (unsigned int i=0; ilagrange_support_points[i]; + this->lagrange_weight = std::sqrt(double(2*k+1)) / prod; } @@ -1307,39 +1161,52 @@ namespace Polynomials } } + + // ------------------ HermiteInterpolation --------------- // HermiteInterpolation::HermiteInterpolation (const unsigned int p) : - Polynomial((p<4) ? 3 : p+1) + Polynomial(0) { + this->coefficients.clear(); + this->in_lagrange_product_form = true; + + this->lagrange_support_points.resize(3); if (p==0) { - this->coefficients[0] = 1.; - this->coefficients[2] = -3.; - this->coefficients[3] = 2.; + this->lagrange_support_points[0] = -0.5; + this->lagrange_support_points[1] = 1.; + this->lagrange_support_points[2] = 1.; + this->lagrange_weight = 2.; } else if (p==1) { - this->coefficients[2] = 3.; - this->coefficients[3] = -2.; + this->lagrange_support_points[0] = 0.; + this->lagrange_support_points[1] = 0.; + this->lagrange_support_points[2] = 1.5; + this->lagrange_weight = -2.; } else if (p==2) { - this->coefficients[1] = 1.; - this->coefficients[2] = -2.; - this->coefficients[3] = 1.; + this->lagrange_support_points[0] = 0.; + this->lagrange_support_points[1] = 1.; + this->lagrange_support_points[2] = 1.; } else if (p==3) { - this->coefficients[2] = -1.; - this->coefficients[3] = 1.; + this->lagrange_support_points[0] = 0.; + this->lagrange_support_points[1] = 0.; + this->lagrange_support_points[2] = 1.; } else { - this->coefficients[4] = 16.; - this->coefficients[3] = -32.; - this->coefficients[2] = 16.; + this->lagrange_support_points.resize(4); + this->lagrange_support_points[0] = 0.; + this->lagrange_support_points[1] = 0.; + this->lagrange_support_points[2] = 1.; + this->lagrange_support_points[3] = 1.; + this->lagrange_weight = 16.; if (p>4) { @@ -1360,6 +1227,7 @@ namespace Polynomials return basis; } + } // ------------------ explicit instantiations --------------- // diff --git a/tests/base/polynomial_legendre_order.cc b/tests/base/polynomial_legendre_order.cc new file mode 100644 index 0000000000..3749d52cf7 --- /dev/null +++ b/tests/base/polynomial_legendre_order.cc @@ -0,0 +1,69 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Similar to polynomial_lagrange_order, but for Legendre interpolation +// This tests the stability of the polynomial evaluation + +#include "../tests.h" +#include +#include +#include + +#include +#include +#include + + +using namespace Polynomials; + + +void check_at_one (const std::vector > &p) +{ + deallog << "Function value of polynomial at right end point: "; + for (unsigned int i=0; i 1e-13*std::sqrt(2*i+1)) + deallog << "Error1 lg y=" << std::log10(std::fabs(y-1.)) + << std::endl; + } + deallog << std::endl; +} + + + +void +check_poly (const unsigned int n) +{ + deallog << "Degree: " << n+1 << std::endl; + std::vector > p = Legendre::generate_complete_basis(n); + check_at_one (p); + deallog << std::endl; +} + + + +int main() +{ + std::ofstream logfile("output"); + deallog << std::setprecision(3); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + check_poly(10); + check_poly(50); +} diff --git a/tests/base/polynomial_legendre_order.output b/tests/base/polynomial_legendre_order.output new file mode 100644 index 0000000000..19cd041f44 --- /dev/null +++ b/tests/base/polynomial_legendre_order.output @@ -0,0 +1,7 @@ + +DEAL::Degree: 11 +DEAL::Function value of polynomial at right end point: ........... +DEAL:: +DEAL::Degree: 51 +DEAL::Function value of polynomial at right end point: ................................................... +DEAL:: diff --git a/tests/serialization/polynomial_legendre.output b/tests/serialization/polynomial_legendre.output index 2681d59721..558ca52379 100644 --- a/tests/serialization/polynomial_legendre.output +++ b/tests/serialization/polynomial_legendre.output @@ -1,4 +1,4 @@ -DEAL::0 0 0 0 4 0 -2.6457513110645907 31.749015732775089 -79.372539331937716 52.915026221291811 0 0 0 1 +DEAL::0 0 0 0 0 0 1 3 0 1.12701665379258298e-01 5.00000000000000000e-01 8.87298334620741702e-01 5.29150262212918179e+01 DEAL::OK