From f41eadf875da6e164db83a59a638a14aaa9fc574 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 9 Feb 2018 18:48:19 +0100 Subject: [PATCH] Add a variant of Hermite polynomials with good condition numbers. --- include/deal.II/base/polynomial.h | 117 +++++++++++++++ source/base/polynomial.cc | 240 ++++++++++++++++++++++++++++++ 2 files changed, 357 insertions(+) diff --git a/include/deal.II/base/polynomial.h b/include/deal.II/base/polynomial.h index ce5fd6772d..6b949900b2 100644 --- a/include/deal.II/base/polynomial.h +++ b/include/deal.II/base/polynomial.h @@ -538,6 +538,7 @@ namespace Polynomials }; + /** * Polynomials for Hermite interpolation condition. * @@ -585,6 +586,122 @@ namespace Polynomials static std::vector > 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: + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + *
 Condition number mass matrix
degreeHermiteInterpolationHermiteLikeInterpolation
n=3105721.40
n=4658015.52
n=51.875e+0418.52
n=66.033e+0419.42
n=109.756e+0527.85
n=159.431e+0640.48
n=252.220e+0868.30
n=352.109e+0998.06
+ * + * 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 + { + public: + /** + * Constructor for the polynomial with index index within the set + * up polynomials of degree @p degree. + */ + HermiteLikeInterpolation (const unsigned int degree, + const unsigned int index); + + /** + * Return the polynomials with index 0 up to degree+1 in + * a space of degree up to degree. Here, degree has to + * be at least 3. + */ + static std::vector > + generate_complete_basis (const unsigned int degree); + }; } diff --git a/source/base/polynomial.cc b/source/base/polynomial.cc index b82971ed1f..dfb805562f 100644 --- a/source/base/polynomial.cc +++ b/source/base/polynomial.cc @@ -18,6 +18,9 @@ #include #include #include +#include +#include + #include #include @@ -1220,6 +1223,9 @@ namespace Polynomials std::vector > HermiteInterpolation::generate_complete_basis (const unsigned int n) { + Assert(n>=3, + ExcNotImplemented("Hermite interpolation makes no sense for " + "degrees less than three")); std::vector > 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 &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(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(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 jacobi_roots(degree-3); + if (degree > 4) + { + LAPACKFullMatrix jacobi_support_points_mat(degree-3, + degree-3); + for (unsigned int k=1; k 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; ilagrange_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; mlagrange_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; mlagrange_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 && indexlagrange_support_points[0] = 0.; + this->lagrange_support_points[1] = 0.; + for (unsigned int m=0, c=2; mlagrange_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; mlagrange_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; mlagrange_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 > + HermiteLikeInterpolation::generate_complete_basis (const unsigned int degree) + { + std::vector > basis (degree + 1); + + for (unsigned int i = 0; i <= degree; ++i) + basis[i] = HermiteLikeInterpolation (degree, i); + + return basis; + } + } // ------------------ explicit instantiations --------------- // -- 2.39.5