]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Switch Legendre polynomials to stable evaluation via root representation 3725/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 3 Jan 2017 17:34:06 +0000 (18:34 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 3 Jan 2017 20:30:20 +0000 (21:30 +0100)
doc/news/changes/minor/20170103MartinKronbichler [new file with mode: 0644]
include/deal.II/base/polynomial.h
source/base/polynomial.cc
tests/base/polynomial_legendre_order.cc [new file with mode: 0644]
tests/base/polynomial_legendre_order.output [new file with mode: 0644]
tests/serialization/polynomial_legendre.output

diff --git a/doc/news/changes/minor/20170103MartinKronbichler b/doc/news/changes/minor/20170103MartinKronbichler
new file mode 100644 (file)
index 0000000..d5453ef
--- /dev/null
@@ -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.
+<br>
+(Martin Kronbichler, 2017/01/03)
index 8e4f22a162680c83f394c639431777e6f4b2f9ad..050b0cf3e7035c5cc7cd30d6c3cfbf94a3815673 100644 (file)
@@ -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 <typename number>
   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 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<number> &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 <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
@@ -390,38 +403,6 @@ namespace Polynomials
     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);
   };
 
   /**
@@ -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 \\
index e921e0dff2ba34a44b9e07edb9a8157db449e6e5..546fea9d0ca7d6068acfd803ee4762256ce72da8 100644 (file)
@@ -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 <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>
@@ -840,177 +841,30 @@ namespace Polynomials
 // ------------------ 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;
   }
 
 
@@ -1307,39 +1161,52 @@ namespace Polynomials
       }
   }
 
+
+
 // ------------------ 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)
           {
@@ -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 (file)
index 0000000..3749d52
--- /dev/null
@@ -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 <iomanip>
+#include <fstream>
+#include <cmath>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/quadrature_lib.h>
+
+
+using namespace Polynomials;
+
+
+void check_at_one (const std::vector<Polynomial<double> > &p)
+{
+  deallog << "Function value of polynomial at right end point: ";
+  for (unsigned int i=0; i<p.size(); ++i)
+    {
+      deallog << '.';
+      const double y = p[i].value(1.);
+      if (std::fabs(y-std::sqrt(2*i+1)) > 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<Polynomial<double> > 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 (file)
index 0000000..19cd041
--- /dev/null
@@ -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::
index 2681d59721e0ddd8b84431642864ea97dd1fb097..558ca5237907da8320c9d3afbdcabdecbe0edeb4 100644 (file)
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.