]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a variant of Hermite polynomials with good condition numbers.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 9 Feb 2018 17:48:19 +0000 (18:48 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 14 Feb 2018 17:28:07 +0000 (18:28 +0100)
include/deal.II/base/polynomial.h
source/base/polynomial.cc

index ce5fd6772d958e6410e35e59821222cf584d6f2d..6b949900b2a8bf952e4b63e806b73f1a6bc658e9 100644 (file)
@@ -538,6 +538,7 @@ namespace Polynomials
   };
 
 
+
   /**
    * Polynomials for Hermite interpolation condition.
    *
@@ -585,6 +586,122 @@ namespace Polynomials
     static std::vector<Polynomial<double> >
     generate_complete_basis (const unsigned int p);
   };
+
+
+
+  /**
+   * Polynomials for a variant of Hermite polynomials with better condition
+   * number in the interpolation than the basis from
+   * HermiteInterpolation. This class is only implemented for degree at least
+   * three, $n\geq 3$.
+   *
+   * In analogy to the actual Hermite polynomials this basis evaluates the
+   * first polynomial $p_0$ to 1 at $x=0$ and has both a zero value and zero
+   * derivative at $x=1$. Likewise, the last polynomial $p_n$ evaluates to 1
+   * at $x=1$ but has zero value and zero derivative at $x=0$. The second
+   * polynomial $p_1$ and the second to last polynomial $p_{n-1}$ represent
+   * the derivative degree of freedom at $x=0$ and $x=1$, respectively. As
+   * such, they are zero at both the end points $x=0, x=1$ and have zero
+   * derivative at the opposite end, $p_1'(1)=0$ and $p_{n-1}'(0)=0$. As
+   * opposed to the original Hermite polynomials, $p_0$ does not have zero
+   * derivative at $x=0$. The additional degree of freedom is used to make
+   * $p_0$ and $p_1$ orthogonal, which for $n=3$ results in a root at
+   * $x=\frac{2}{7}$ for $p_0$ and at $x=\frac{5}{7}$ for $p_n$,
+   * respectively. Furthermore, the extension of these polynomials to higher
+   * degrees $n>3$ is constructed by adding additional nodes inside the unit
+   * interval, again ensuring better conditioning. The nodes are computed as
+   * the roots of the Jacobi polynomials for $\alpha=\beta=2$ which are
+   * orthogonal against the generating function $x^2(1-x)^2$ with the Hermite
+   * property. Then, these polynomials are constructed in the usual way as
+   * Lagrange polynomials with double roots at $x=0$ and $x=1$. For example at
+   * $n=4$, all of $p_0, p_1, p_3, p_4$ get an additional root at $x=0.5$
+   * through the factor $(x-0.5)$.
+
+   * These two relaxations improve the condition number of the mass matrix
+   * (i.e., interpolation) significantly, as can be seen from the following
+   * table:
+   *
+   * <table align="center" border="1">
+   *   <tr>
+   *    <th>&nbsp;</th>
+   *    <th colspan="2">Condition number mass matrix</th>
+   *   </tr>
+   *   <tr>
+   *    <th>degree</th>
+   *    <th>HermiteInterpolation</th>
+   *    <th>HermiteLikeInterpolation</th>
+   *   </tr>
+   *   <tr>
+   *    <th>n=3</th>
+   *    <th>1057</th>
+   *    <th>21.40</th>
+   *   </tr>
+   *   <tr>
+   *    <th>n=4</th>
+   *    <th>6580</th>
+   *    <th>15.52</th>
+   *   </tr>
+   *   <tr>
+   *    <th>n=5</th>
+   *    <th>1.875e+04</th>
+   *    <th>18.52</th>
+   *   </tr>
+   *   <tr>
+   *    <th>n=6</th>
+   *    <th>6.033e+04</th>
+   *    <th>19.42</th>
+   *   </tr>
+   *   <tr>
+   *    <th>n=10</th>
+   *    <th>9.756e+05</th>
+   *    <th>27.85</th>
+   *   </tr>
+   *   <tr>
+   *    <th>n=15</th>
+   *    <th>9.431e+06</th>
+   *    <th>40.48</th>
+   *   </tr>
+   *   <tr>
+   *    <th>n=25</th>
+   *    <th>2.220e+08</th>
+   *    <th>68.30</th>
+   *   </tr>
+   *   <tr>
+   *    <th>n=35</th>
+   *    <th>2.109e+09</th>
+   *    <th>98.06</th>
+   *   </tr>
+   * </table>
+   *
+   * This polynomial inherits the advantageous property of Hermite polynomials
+   * where only two functions have value and/or derivative nonzero on a face
+   * but gives better condition numbers of interpolation, which improves the
+   * performance of some iterative schemes like conjugate gradients with
+   * point-Jacobi.
+   *
+   * @note This class requires LAPACK support.
+   *
+   * @author Martin Kronbichler
+   * @date 2018
+   */
+  class HermiteLikeInterpolation : public Polynomial<double>
+  {
+  public:
+    /**
+     * Constructor for the polynomial with index <tt>index</tt> within the set
+     * up polynomials of degree @p degree.
+     */
+    HermiteLikeInterpolation (const unsigned int degree,
+                              const unsigned int index);
+
+    /**
+     * Return the polynomials with index <tt>0</tt> up to <tt>degree+1</tt> in
+     * a space of degree up to <tt>degree</tt>. Here, <tt>degree</tt> has to
+     * be at least 3.
+     */
+    static std::vector<Polynomial<double> >
+    generate_complete_basis (const unsigned int degree);
+  };
 }
 
 
index b82971ed1f034ad0007db9f36ec30f78bf4545af..dfb805562f9bab5c67f88a4a4d9b9c96a360ea08 100644 (file)
@@ -18,6 +18,9 @@
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/thread_management.h>
 #include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/lapack_full_matrix.h>
+
 
 #include <cmath>
 #include <algorithm>
@@ -1220,6 +1223,9 @@ namespace Polynomials
   std::vector<Polynomial<double> >
   HermiteInterpolation::generate_complete_basis (const unsigned int n)
   {
+    Assert(n>=3,
+           ExcNotImplemented("Hermite interpolation makes no sense for "
+                             "degrees less than three"));
     std::vector<Polynomial<double> > basis (n + 1);
 
     for (unsigned int i = 0; i <= n; ++i)
@@ -1228,6 +1234,240 @@ namespace Polynomials
     return basis;
   }
 
+
+// ------------------ HermiteLikeInterpolation --------------- //
+  namespace
+  {
+    // Finds the zero position x_star such that the mass matrix entry (0,1)
+    // with the Hermite polynomials evaluates to zero. The function has
+    // originally been derived by a secant method for the integral entry
+    // l_0(x) * l_1(x) but we only need to do one iteration because the zero
+    // x_star is linear in the integral value.
+    double find_support_point_x_star (const Vector<double> &jacobi_roots)
+    {
+      // Initial guess for the support point position values: The zero turns
+      // out to be between zero and the first root of the Jacobi polynomial,
+      // but the algorithm is agnostic about that, so simply choose two points
+      // that are sufficiently far apart.
+      double guess_left = 0;
+      double guess_right = 0.5;
+      const unsigned int degree = jacobi_roots.size() + 3;
+
+      // Compute two integrals of the product of l_0(x) * l_1(x)
+      // l_0(x) = (x-y)*(x-jacobi_roots(0))*...*(x-jacobi_roos(degree-4))*(x-1)*(x-1)
+      // l_1(x) = (x-0)*(x-jacobi_roots(0))*...*(x-jacobi_roots(degree-4))*(x-1)*(x-1)
+      // where y is either guess_left or guess_right for the two integrals.
+      // Note that the polynomials are not yet normalized here, which is not
+      // necessary because we are only looking for the x_star where the matrix
+      // entry is zero, for which the constants do not matter.
+      QGauss<1> gauss(degree + 1);
+      double integral_left = 0, integral_right = 0;
+      for (unsigned int q=0; q<gauss.size(); ++q)
+        {
+          const double x = gauss.point(q)[0];
+          double poly_val_common = x;
+          for (unsigned int j=0; j<degree-3; ++j)
+            poly_val_common *= Utilities::fixed_power<2>(x-jacobi_roots(j));
+          poly_val_common *= Utilities::fixed_power<4>(x - 1.);
+          integral_left += gauss.weight(q)*(poly_val_common*(x - guess_left));
+          integral_right += gauss.weight(q)*(poly_val_common*(x - guess_right));
+        }
+
+      // compute guess by secant method. Due to linearity in the root x_star,
+      // this is the correct position after this single step
+      return guess_right - (guess_right-guess_left)/(integral_right-integral_left)*integral_right;
+    }
+  }
+
+
+
+  HermiteLikeInterpolation::HermiteLikeInterpolation (const unsigned int degree,
+                                                      const unsigned int index)
+    :
+    Polynomial<double>(0)
+  {
+    Assert(degree>=3,
+           ExcNotImplemented("Hermite interpolation makes no sense for "
+                             "degrees less than three"));
+    AssertIndexRange(index, degree+1);
+
+    this->coefficients.clear();
+    this->in_lagrange_product_form = true;
+
+    this->lagrange_support_points.resize(degree);
+
+    // 4 Polynomials with degree 3
+    // entries (1,0) and (3,2) of the mass matrix will be equal to 0
+    //
+    //     | x  0  x  x |
+    //     | 0  x  x  x |
+    // M = | x  x  x  0 |
+    //     | x  x  0  x |
+    //
+    if (degree==3)
+      {
+        if (index==0)
+          {
+            this->lagrange_support_points[0] = 2./7.;
+            this->lagrange_support_points[1] = 1.;
+            this->lagrange_support_points[2] = 1.;
+            this->lagrange_weight = -3.5;
+          }
+        else if (index==1)
+          {
+            this->lagrange_support_points[0] = 0.;
+            this->lagrange_support_points[1] = 1.;
+            this->lagrange_support_points[2] = 1.;
+            this->lagrange_weight = 6.75;
+          }
+        else if (index==2)
+          {
+            this->lagrange_support_points[0] = 0.;
+            this->lagrange_support_points[1] = 0.;
+            this->lagrange_support_points[2] = 1.;
+            this->lagrange_weight = -6.75;
+          }
+        else if (index==3)
+          {
+            this->lagrange_support_points[0] = 0.;
+            this->lagrange_support_points[1] = 0.;
+            this->lagrange_support_points[2] = 5./7.;
+            this->lagrange_weight = 3.5;
+          }
+      }
+
+    // Higher order Polynomials degree>=4: the entries (1,0) and
+    // (degree,degree-1) of the mass matrix will be equal to 0
+    //
+    //     | x  0  x  x         x  x  x |
+    //     | 0  x  x  x  . . .  x  x  x |
+    //     | x  x  x  x         x  x  x |
+    //     | x  x  x  x         x  x  x |
+    //     |     .       .         .    |
+    // M = |     .         .       .    |
+    //     |     .           .     .    |
+    //     | x  x  x  x         x  x  x |
+    //     | x  x  x  x  . . .  x  x  0 |
+    //     | x  x  x  x         x  0  x |
+    //
+    if (degree >= 4)
+      {
+        // We find the inner points as the zeros of the Jacobi polynomials
+        // with alpha = beta = 2 which is the polynomial with the kernel
+        // (1-x)^2 (1+x)^2, the two polynomials achieving zero value and zero
+        // derivative at the boundary. The zeros of the Jacobi polynomials are
+        // given by the eigenvalues to a symmetric tridiagonal matrix with the
+        // entries given below. For degree 4 the eigenvalue is zero, so bypass
+        // the LAPACK logic in that case.
+
+        Vector<double> jacobi_roots(degree-3);
+        if (degree > 4)
+          {
+            LAPACKFullMatrix<double> jacobi_support_points_mat(degree-3,
+                                                               degree-3);
+            for (unsigned int k=1; k<degree-3; k++)
+              {
+                jacobi_support_points_mat(k-1,k) =
+                  std::sqrt(4.*k*(k+2.)*(k+2.)*(k+4.)/((2.*k+3.)*(2.*k+4.)*(2.*k+4.)*(2.*k+5.)));
+                jacobi_support_points_mat(k,k-1) = jacobi_support_points_mat(k-1,k);
+              }
+
+            // calculate the eigenvalues = zero points of the Jacobi polynomials
+            FullMatrix<double> eigenvectors(degree-3,degree-3);
+            jacobi_support_points_mat.compute_eigenvalues_symmetric(-1., 1., 1.e-20,
+                                                                    jacobi_roots,
+                                                                    eigenvectors);
+            AssertDimension(jacobi_roots.size(), degree-3);
+
+            // Note that this algorithm computes the zeros of the Jacobi
+            // polynomial for the interval [-1,1], so we must scale the
+            // eigenvalues to the interval [0,1] before using them
+            for (unsigned int i=0; i<degree-3; ++i)
+              jacobi_roots(i) = 0.5*jacobi_roots(i)+0.5;
+          }
+        else
+          // only a single zero at x=0.5 for degree==4
+          jacobi_roots(0) = 0.5;
+
+        // iteration from variable support point N with secant method
+        // initial values
+
+        this->lagrange_support_points.resize(degree);
+        if (index==0)
+          {
+            const double auxiliary_zero = find_support_point_x_star(jacobi_roots);
+            this->lagrange_support_points[0] = auxiliary_zero;
+            for (unsigned int m=0; m<degree-3; m++)
+              this->lagrange_support_points[m+1] = jacobi_roots(m);
+            this->lagrange_support_points[degree-2] = 1.;
+            this->lagrange_support_points[degree-1] = 1.;
+
+            // ensure that the polynomial evaluates to one at x=0
+            this->lagrange_weight = 1./this->value(0);
+          }
+        else if (index==1)
+          {
+            this->lagrange_support_points[0] = 0.;
+            for (unsigned int m=0; m<degree-3; m++)
+              this->lagrange_support_points[m+1] = jacobi_roots(m);
+            this->lagrange_support_points[degree-2] = 1.;
+            this->lagrange_support_points[degree-1] = 1.;
+
+            // scale close to approximate maximum
+            this->lagrange_weight = 1./this->value(0.4*jacobi_roots(0));
+          }
+        else if (index>=2 && index<degree-1)
+          {
+            this->lagrange_support_points[0] = 0.;
+            this->lagrange_support_points[1] = 0.;
+            for (unsigned int m=0, c=2; m<degree-3; m++)
+              if (m+2 != index)
+                this->lagrange_support_points[c++] = jacobi_roots(m);
+            this->lagrange_support_points[degree-2] = 1.;
+            this->lagrange_support_points[degree-1] = 1.;
+
+            // ensure that the polynomial evaluates to one at the respective
+            // nodal point
+            this->lagrange_weight = 1./this->value(jacobi_roots(index-2));
+          }
+        else if (index==degree-1)
+          {
+            this->lagrange_support_points[0] = 0.;
+            this->lagrange_support_points[1] = 0.;
+            for (unsigned int m=0; m<degree-3; m++)
+              this->lagrange_support_points[m+2] = jacobi_roots(m);
+            this->lagrange_support_points[degree-1] = 1.;
+
+            this->lagrange_weight = 1./this->value(1.0-0.4*jacobi_roots(0));
+          }
+        else if (index==degree)
+          {
+            const double auxiliary_zero = find_support_point_x_star(jacobi_roots);
+            this->lagrange_support_points[0] = 0.;
+            this->lagrange_support_points[1] = 0.;
+            for (unsigned int m=0; m<degree-3; m++)
+              this->lagrange_support_points[m+2] = jacobi_roots(m);
+            this->lagrange_support_points[degree-1] = 1.-auxiliary_zero;
+
+            // ensure that the polynomial evaluates to one at x=1
+            this->lagrange_weight = 1./this->value(1.);
+          }
+      }
+  }
+
+
+
+  std::vector<Polynomial<double> >
+  HermiteLikeInterpolation::generate_complete_basis (const unsigned int degree)
+  {
+    std::vector<Polynomial<double> > basis (degree + 1);
+
+    for (unsigned int i = 0; i <= degree; ++i)
+      basis[i] = HermiteLikeInterpolation (degree, i);
+
+    return basis;
+  }
+
 }
 
 // ------------------ explicit instantiations --------------- //

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.