// ---------------------------------------------------------------------
//
-// 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.
//
/**
* 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 <typename number>
class Polynomial : public Subscriptor
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
* 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;
* 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<number> &values) const;
* 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,
/**
* Legendre polynomials of arbitrary degree. Constructing a Legendre
- * polynomial of degree <tt>p</tt>, the coefficients will be computed by the
- * three-term recursion formula.
+ * polynomial of degree <tt>p</tt>, 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
static
std::vector<Polynomial<double> >
generate_complete_basis (const unsigned int degree);
-
- private:
- /**
- * Coefficients for the interval $[0,1]$.
- */
- static std::vector<std_cxx11::shared_ptr<const std::vector<double> > > 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<std_cxx11::shared_ptr<const std::vector<double> > > 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 <tt>scale</tt> and <tt>shift</tt> functions
- * of <tt>Polynomial</tt>, 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<double> &
- get_coefficients (const unsigned int k);
};
/**
* 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 \\
// ---------------------------------------------------------------------
//
-// 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.
//
#include <deal.II/base/point.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/thread_management.h>
+#include <deal.II/base/quadrature_lib.h>
#include <cmath>
#include <algorithm>
// ------------------ class Legendre --------------- //
-// Reserve space for polynomials up to degree 19. Should be sufficient
-// for the start.
- std::vector<std_cxx11::shared_ptr<const std::vector<double> > >
- Legendre::recursive_coefficients(20);
- std::vector<std_cxx11::shared_ptr<const std::vector<double> > >
- Legendre::shifted_coefficients(20);
-
Legendre::Legendre (const unsigned int k)
:
- Polynomial<double> (get_coefficients(k))
- {}
-
-
-
- void
- Legendre::compute_coefficients (const unsigned int k_)
+ Polynomial<double> (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<const std::vector<double> >())))
- // 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<double> *c0 = new std::vector<double>(1);
- (*c0)[0] = 1.;
-
- std::vector<double> *c1 = new std::vector<double>(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<const std::vector<double> >(c0);
- recursive_coefficients[1] =
- std_cxx11::shared_ptr<const std::vector<double> >(c1);
-
- // Compute polynomials
- // orthogonal on [0,1]
- c0 = new std::vector<double>(*c0);
- c1 = new std::vector<double>(*c1);
-
- Polynomial<double>::shift<SHIFT_TYPE> (*c0, -1.);
- Polynomial<double>::scale(*c0, 2.);
- Polynomial<double>::shift<SHIFT_TYPE> (*c1, -1.);
- Polynomial<double>::scale(*c1, 2.);
- Polynomial<double>::multiply(*c1, std::sqrt(3.));
- shifted_coefficients[0]=std_cxx11::shared_ptr<const std::vector<double> >(c0);
- shifted_coefficients[1]=std_cxx11::shared_ptr<const std::vector<double> >(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<double> *ck = new std::vector<double>(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<const std::vector<double> >(ck);
- // and compute the
- // coefficients for [0,1]
- ck = new std::vector<double>(*ck);
- Polynomial<double>::shift<SHIFT_TYPE> (*ck, -1.);
- Polynomial<double>::scale(*ck, 2.);
- Polynomial<double>::multiply(*ck, std::sqrt(2.*k+1.));
- shifted_coefficients[k] =
- std_cxx11::shared_ptr<const std::vector<double> >(ck);
- };
- };
- }
-
-
-
- const std::vector<double> &
- 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; i<k; ++i)
+ this->lagrange_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; i<k; ++i)
+ prod *= this->lagrange_support_points[i];
+ this->lagrange_weight = std::sqrt(double(2*k+1)) / prod;
}
}
}
+
+
// ------------------ HermiteInterpolation --------------- //
HermiteInterpolation::HermiteInterpolation (const unsigned int p)
:
- Polynomial<double>((p<4) ? 3 : p+1)
+ Polynomial<double>(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)
{
return basis;
}
+
}
// ------------------ explicit instantiations --------------- //
projects / dealii.git / commitdiff