Reference documentation for deal.II version 9.4.0
polynomial.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2000 - 2021 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_polynomial_h
17#define dealii_polynomial_h
18
19
20
21#include <deal.II/base/config.h>
22
24#include <deal.II/base/point.h>
26
27#include <array>
28#include <memory>
29#include <vector>
30
32
42namespace Polynomials
43{
63 template <typename number>
64 class Polynomial : public Subscriptor
65 {
66 public:
75 Polynomial(const std::vector<number> &coefficients);
76
80 Polynomial(const unsigned int n);
81
90 const unsigned int evaluation_point);
91
96
106 number
107 value(const number x) const;
108
119 void
120 value(const number x, std::vector<number> &values) const;
121
140 template <typename Number2>
141 void
142 value(const Number2 x,
143 const unsigned int n_derivatives,
144 Number2 * values) const;
145
158 template <std::size_t n_entries, typename Number2>
159 void
160 values_of_array(const std::array<Number2, n_entries> &points,
161 const unsigned int n_derivatives,
162 std::array<Number2, n_entries> * values) const;
163
169 unsigned int
170 degree() const;
171
179 void
180 scale(const number factor);
181
197 template <typename number2>
198 void
199 shift(const number2 offset);
200
205 derivative() const;
206
212 primitive() const;
213
218 operator*=(const double s);
219
225
231
237
241 bool
242 operator==(const Polynomial<number> &p) const;
243
247 void
248 print(std::ostream &out) const;
249
255 template <class Archive>
256 void
257 serialize(Archive &ar, const unsigned int version);
258
262 virtual std::size_t
263 memory_consumption() const;
264
265 protected:
269 static void
270 scale(std::vector<number> &coefficients, const number factor);
271
275 template <typename number2>
276 static void
277 shift(std::vector<number> &coefficients, const number2 shift);
278
282 static void
283 multiply(std::vector<number> &coefficients, const number factor);
284
290 void
292
301 std::vector<number> coefficients;
302
308
313 std::vector<number> lagrange_support_points;
314
320 };
321
322
327 template <typename number>
328 class Monomial : public Polynomial<number>
329 {
330 public:
335 Monomial(const unsigned int n, const double coefficient = 1.);
336
343 static std::vector<Polynomial<number>>
344 generate_complete_basis(const unsigned int degree);
345
346 private:
350 static std::vector<number>
351 make_vector(unsigned int n, const double coefficient);
352 };
353
354
370 class LagrangeEquidistant : public Polynomial<double>
371 {
372 public:
378 LagrangeEquidistant(const unsigned int n, const unsigned int support_point);
379
388 static std::vector<Polynomial<double>>
389 generate_complete_basis(const unsigned int degree);
390
391 private:
396 static void
397 compute_coefficients(const unsigned int n,
398 const unsigned int support_point,
399 std::vector<double> &a);
400 };
401
402
403
410 std::vector<Polynomial<double>>
411 generate_complete_Lagrange_basis(const std::vector<Point<1>> &points);
412
413
414
428 class Legendre : public Polynomial<double>
429 {
430 public:
434 Legendre(const unsigned int p);
435
442 static std::vector<Polynomial<double>>
443 generate_complete_basis(const unsigned int degree);
444 };
445
465 class Lobatto : public Polynomial<double>
466 {
467 public:
472 Lobatto(const unsigned int p = 0);
473
478 static std::vector<Polynomial<double>>
479 generate_complete_basis(const unsigned int p);
480
481 private:
485 std::vector<double>
486 compute_coefficients(const unsigned int p);
487 };
488
489
490
528 class Hierarchical : public Polynomial<double>
529 {
530 public:
535 Hierarchical(const unsigned int p);
536
547 static std::vector<Polynomial<double>>
548 generate_complete_basis(const unsigned int degree);
549
550 private:
554 static void
555 compute_coefficients(const unsigned int p);
556
561 static const std::vector<double> &
562 get_coefficients(const unsigned int p);
563
572 static std::vector<std::unique_ptr<const std::vector<double>>>
574 };
575
576
577
605 class HermiteInterpolation : public Polynomial<double>
606 {
607 public:
612 HermiteInterpolation(const unsigned int p);
613
619 static std::vector<Polynomial<double>>
620 generate_complete_basis(const unsigned int p);
621 };
622
623
624
727 {
728 public:
733 HermiteLikeInterpolation(const unsigned int degree,
734 const unsigned int index);
735
740 static std::vector<Polynomial<double>>
741 generate_complete_basis(const unsigned int degree);
742 };
743
744
745
746 /*
747 * Evaluate a Jacobi polynomial @f$P_n^{\alpha, \beta}(x) @f$ specified by the
748 * parameters @p alpha, @p beta, @p n, where @p n is the degree of the
749 * Jacobi polynomial.
750 *
751 * @note The Jacobi polynomials are not orthonormal and are defined on the
752 * unit interval @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often
753 * used in literature. @p x is the point of evaluation.
754 */
755 template <typename Number>
756 Number
757 jacobi_polynomial_value(const unsigned int degree,
758 const int alpha,
759 const int beta,
760 const Number x);
761
762
775 template <typename Number>
776 std::vector<Number>
777 jacobi_polynomial_roots(const unsigned int degree,
778 const int alpha,
779 const int beta);
780} // namespace Polynomials
781
782
785/* -------------------------- inline functions --------------------- */
786
787namespace Polynomials
788{
789 template <typename number>
791 : in_lagrange_product_form(false)
792 , lagrange_weight(1.)
793 {}
794
795
796
797 template <typename number>
798 inline unsigned int
800 {
801 if (in_lagrange_product_form == true)
802 {
803 return lagrange_support_points.size();
804 }
805 else
806 {
807 Assert(coefficients.size() > 0, ExcEmptyObject());
808 return coefficients.size() - 1;
809 }
810 }
811
812
813
814 template <typename number>
815 inline number
816 Polynomial<number>::value(const number x) const
817 {
818 if (in_lagrange_product_form == false)
819 {
820 Assert(coefficients.size() > 0, ExcEmptyObject());
821
822 // Horner scheme
823 const unsigned int m = coefficients.size();
824 number value = coefficients.back();
825 for (int k = m - 2; k >= 0; --k)
826 value = value * x + coefficients[k];
827 return value;
828 }
829 else
830 {
831 // direct evaluation of Lagrange polynomial
832 const unsigned int m = lagrange_support_points.size();
833 number value = 1.;
834 for (unsigned int j = 0; j < m; ++j)
835 value *= x - lagrange_support_points[j];
836 value *= lagrange_weight;
837 return value;
838 }
839 }
840
841
842
843 template <typename number>
844 template <typename Number2>
845 inline void
847 const unsigned int n_derivatives,
848 Number2 * values) const
849 {
850 values_of_array(std::array<Number2, 1ul>{{x}},
851 n_derivatives,
852 reinterpret_cast<std::array<Number2, 1ul> *>(values));
853 }
854
855
856
857 template <typename number>
858 template <std::size_t n_entries, typename Number2>
859 inline void
861 const std::array<Number2, n_entries> &x,
862 const unsigned int n_derivatives,
863 std::array<Number2, n_entries> * values) const
864 {
865 // evaluate Lagrange polynomial and derivatives
866 if (in_lagrange_product_form == true)
867 {
868 // to compute the value and all derivatives of a polynomial of the
869 // form (x-x_1)*(x-x_2)*...*(x-x_n), expand the derivatives like
870 // automatic differentiation does.
871 const unsigned int n_supp = lagrange_support_points.size();
872 const number weight = lagrange_weight;
873 switch (n_derivatives)
874 {
875 default:
876 for (unsigned int e = 0; e < n_entries; ++e)
877 values[0][e] = 1.;
878 for (unsigned int k = 1; k <= n_derivatives; ++k)
879 for (unsigned int e = 0; e < n_entries; ++e)
880 values[k][e] = 0.;
881 for (unsigned int i = 0; i < n_supp; ++i)
882 {
883 std::array<Number2, n_entries> v = x;
884 for (unsigned int e = 0; e < n_entries; ++e)
885 v[e] -= lagrange_support_points[i];
886
887 // multiply by (x-x_i) and compute action on all derivatives,
888 // too (inspired from automatic differentiation: implement the
889 // product rule for the old value and the new variable 'v',
890 // i.e., expand value v and derivative one). since we reuse a
891 // value from the next lower derivative from the steps before,
892 // need to start from the highest derivative
893 for (unsigned int k = n_derivatives; k > 0; --k)
894 for (unsigned int e = 0; e < n_entries; ++e)
895 values[k][e] = (values[k][e] * v[e] + values[k - 1][e]);
896 for (unsigned int e = 0; e < n_entries; ++e)
897 values[0][e] *= v[e];
898 }
899 // finally, multiply by the weight in the Lagrange
900 // denominator. Could be done instead of setting values[0] = 1
901 // above, but that gives different accumulation of round-off
902 // errors (multiplication is not associative) compared to when we
903 // computed the weight, and hence a basis function might not be
904 // exactly one at the center point, which is nice to have. We also
905 // multiply derivatives by k! to transform the product p_n =
906 // p^(n)(x)/k! into the actual form of the derivative
907 {
908 number k_factorial = 1;
909 for (unsigned int k = 0; k <= n_derivatives; ++k)
910 {
911 for (unsigned int e = 0; e < n_entries; ++e)
912 values[k][e] *= k_factorial * weight;
913 k_factorial *= static_cast<number>(k + 1);
914 }
915 }
916 break;
917
918 // manually implement case 0 (values only), case 1 (value + first
919 // derivative), and case 2 (up to second derivative) since they
920 // might be called often. then, we can unroll the inner loop and
921 // keep the temporary results as local variables to help the
922 // compiler with the pointer aliasing analysis.
923 case 0:
924 {
925 std::array<Number2, n_entries> value;
926 for (unsigned int e = 0; e < n_entries; ++e)
927 value[e] = 1.;
928 for (unsigned int i = 0; i < n_supp; ++i)
929 for (unsigned int e = 0; e < n_entries; ++e)
930 value[e] *= (x[e] - lagrange_support_points[i]);
931
932 for (unsigned int e = 0; e < n_entries; ++e)
933 values[0][e] = weight * value[e];
934 break;
935 }
936
937 case 1:
938 {
939 std::array<Number2, n_entries> value, derivative = {};
940 for (unsigned int e = 0; e < n_entries; ++e)
941 value[e] = 1.;
942 for (unsigned int i = 0; i < n_supp; ++i)
943 for (unsigned int e = 0; e < n_entries; ++e)
944 {
945 const Number2 v = x[e] - lagrange_support_points[i];
946 derivative[e] = derivative[e] * v + value[e];
947 value[e] *= v;
948 }
949
950 for (unsigned int e = 0; e < n_entries; ++e)
951 {
952 values[0][e] = weight * value[e];
953 values[1][e] = weight * derivative[e];
954 }
955 break;
956 }
957
958 case 2:
959 {
960 std::array<Number2, n_entries> value, derivative = {},
961 second = {};
962 for (unsigned int e = 0; e < n_entries; ++e)
963 value[e] = 1.;
964 for (unsigned int i = 0; i < n_supp; ++i)
965 for (unsigned int e = 0; e < n_entries; ++e)
966 {
967 const Number2 v = x[e] - lagrange_support_points[i];
968 second[e] = second[e] * v + derivative[e];
969 derivative[e] = derivative[e] * v + value[e];
970 value[e] *= v;
971 }
972
973 for (unsigned int e = 0; e < n_entries; ++e)
974 {
975 values[0][e] = weight * value[e];
976 values[1][e] = weight * derivative[e];
977 values[2][e] = static_cast<number>(2) * weight * second[e];
978 }
979 break;
980 }
981 }
982 return;
983 }
984
985 Assert(coefficients.size() > 0, ExcEmptyObject());
986
987 // if derivatives are needed, then do it properly by the full
988 // Horner scheme
989 const unsigned int m = coefficients.size();
990 std::vector<std::array<Number2, n_entries>> a(coefficients.size());
991 for (unsigned int i = 0; i < coefficients.size(); ++i)
992 for (unsigned int e = 0; e < n_entries; ++e)
993 a[i][e] = coefficients[i];
994
995 unsigned int j_factorial = 1;
996
997 // loop over all requested derivatives. note that derivatives @p{j>m} are
998 // necessarily zero, as they differentiate the polynomial more often than
999 // the highest power is
1000 const unsigned int min_valuessize_m = std::min(n_derivatives + 1, m);
1001 for (unsigned int j = 0; j < min_valuessize_m; ++j)
1002 {
1003 for (int k = m - 2; k >= static_cast<int>(j); --k)
1004 for (unsigned int e = 0; e < n_entries; ++e)
1005 a[k][e] += x[e] * a[k + 1][e];
1006 for (unsigned int e = 0; e < n_entries; ++e)
1007 values[j][e] = static_cast<number>(j_factorial) * a[j][e];
1008
1009 j_factorial *= j + 1;
1010 }
1011
1012 // fill higher derivatives by zero
1013 for (unsigned int j = min_valuessize_m; j <= n_derivatives; ++j)
1014 for (unsigned int e = 0; e < n_entries; ++e)
1015 values[j][e] = 0.;
1016 }
1017
1018
1019
1020 template <typename number>
1021 template <class Archive>
1022 inline void
1023 Polynomial<number>::serialize(Archive &ar, const unsigned int)
1024 {
1025 // forward to serialization function in the base class.
1026 ar &static_cast<Subscriptor &>(*this);
1027 ar &coefficients;
1028 ar &in_lagrange_product_form;
1029 ar &lagrange_support_points;
1030 ar &lagrange_weight;
1031 }
1032
1033
1034
1035 template <typename Number>
1036 Number
1037 jacobi_polynomial_value(const unsigned int degree,
1038 const int alpha,
1039 const int beta,
1040 const Number x)
1041 {
1042 Assert(alpha >= 0 && beta >= 0,
1043 ExcNotImplemented("Negative alpha/beta coefficients not supported"));
1044 // the Jacobi polynomial is evaluated using a recursion formula.
1045 Number p0, p1;
1046
1047 // The recursion formula is defined for the interval [-1, 1], so rescale
1048 // to that interval here
1049 const Number xeval = Number(-1) + 2. * x;
1050
1051 // initial values P_0(x), P_1(x):
1052 p0 = 1.0;
1053 if (degree == 0)
1054 return p0;
1055 p1 = ((alpha + beta + 2) * xeval + (alpha - beta)) / 2;
1056 if (degree == 1)
1057 return p1;
1058
1059 for (unsigned int i = 1; i < degree; ++i)
1060 {
1061 const Number v = 2 * i + (alpha + beta);
1062 const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
1063 const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
1064 const Number a3 = v * (v + 1) * (v + 2);
1065 const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
1066
1067 const Number pn = ((a2 + a3 * xeval) * p1 - a4 * p0) / a1;
1068 p0 = p1;
1069 p1 = pn;
1070 }
1071 return p1;
1072 }
1073
1074
1075
1076 template <typename Number>
1077 std::vector<Number>
1078 jacobi_polynomial_roots(const unsigned int degree,
1079 const int alpha,
1080 const int beta)
1081 {
1082 std::vector<Number> x(degree, 0.5);
1083
1084 // compute zeros with a Newton algorithm.
1085
1086 // Set tolerance. For long double we might not always get the additional
1087 // precision in a run time environment (e.g. with valgrind), so we must
1088 // limit the tolerance to double. Since we do a Newton iteration, doing
1089 // one more iteration after the residual has indicated convergence will be
1090 // enough for all number types due to the quadratic convergence of
1091 // Newton's method
1092
1093 const Number tolerance =
1094 4 * std::max(static_cast<Number>(std::numeric_limits<double>::epsilon()),
1096
1097 // The following implementation follows closely the one given in the
1098 // appendix of the book by Karniadakis and Sherwin: Spectral/hp element
1099 // methods for computational fluid dynamics (Oxford University Press,
1100 // 2005)
1101
1102 // If symmetric, we only need to compute the half of points
1103 const unsigned int n_points = (alpha == beta ? degree / 2 : degree);
1104 for (unsigned int k = 0; k < n_points; ++k)
1105 {
1106 // we take the zeros of the Chebyshev polynomial (alpha=beta=-0.5) as
1107 // initial values, corrected by the initial value
1108 Number r = 0.5 - 0.5 * std::cos(static_cast<Number>(2 * k + 1) /
1109 (2 * degree) * numbers::PI);
1110 if (k > 0)
1111 r = (r + x[k - 1]) / 2;
1112
1113 unsigned int converged = numbers::invalid_unsigned_int;
1114 for (unsigned int it = 1; it < 1000; ++it)
1115 {
1116 Number s = 0.;
1117 for (unsigned int i = 0; i < k; ++i)
1118 s += 1. / (r - x[i]);
1119
1120 // derivative of P_n^{alpha,beta}, rescaled to [0, 1]
1121 const Number J_x =
1122 (alpha + beta + degree + 1) *
1123 jacobi_polynomial_value(degree - 1, alpha + 1, beta + 1, r);
1124
1125 // value of P_n^{alpha,beta}
1126 const Number f = jacobi_polynomial_value(degree, alpha, beta, r);
1127 const Number delta = f / (f * s - J_x);
1128 r += delta;
1129 if (converged == numbers::invalid_unsigned_int &&
1130 std::abs(delta) < tolerance)
1131 converged = it;
1132
1133 // do one more iteration to ensure accuracy also for tighter
1134 // types than double (e.g. long double)
1135 if (it == converged + 1)
1136 break;
1137 }
1138
1140 ExcMessage("Newton iteration for zero of Jacobi polynomial "
1141 "did not converge."));
1142
1143 x[k] = r;
1144 }
1145
1146 // in case we assumed symmetry, fill up the missing values
1147 for (unsigned int k = n_points; k < degree; ++k)
1148 x[k] = 1.0 - x[degree - k - 1];
1149
1150 return x;
1151 }
1152
1153} // namespace Polynomials
1155
1156#endif
HermiteInterpolation(const unsigned int p)
Definition: polynomial.cc:1057
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
Definition: polynomial.cc:1109
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:1405
HermiteLikeInterpolation(const unsigned int degree, const unsigned int index)
Definition: polynomial.cc:1175
static std::vector< std::unique_ptr< const std::vector< double > > > recursive_coefficients
Definition: polynomial.h:573
Hierarchical(const unsigned int p)
Definition: polynomial.cc:872
static void compute_coefficients(const unsigned int p)
Definition: polynomial.cc:879
static const std::vector< double > & get_coefficients(const unsigned int p)
Definition: polynomial.cc:1013
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:1029
static void compute_coefficients(const unsigned int n, const unsigned int support_point, std::vector< double > &a)
Definition: polynomial.cc:621
LagrangeEquidistant(const unsigned int n, const unsigned int support_point)
Definition: polynomial.cc:597
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:679
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:745
Legendre(const unsigned int p)
Definition: polynomial.cc:718
std::vector< double > compute_coefficients(const unsigned int p)
Definition: polynomial.cc:764
Lobatto(const unsigned int p=0)
Definition: polynomial.cc:759
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
Definition: polynomial.cc:850
static std::vector< Polynomial< number > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:566
Monomial(const unsigned int n, const double coefficient=1.)
Definition: polynomial.cc:558
static std::vector< number > make_vector(unsigned int n, const double coefficient)
Definition: polynomial.cc:548
number value(const number x) const
Definition: polynomial.h:816
bool operator==(const Polynomial< number > &p) const
Definition: polynomial.cc:347
std::vector< number > coefficients
Definition: polynomial.h:301
Polynomial< number > primitive() const
Definition: polynomial.cc:488
Polynomial< number > & operator+=(const Polynomial< number > &p)
Definition: polynomial.cc:269
void values_of_array(const std::array< Number2, n_entries > &points, const unsigned int n_derivatives, std::array< Number2, n_entries > *values) const
Definition: polynomial.h:860
Polynomial< number > derivative() const
Definition: polynomial.cc:459
void transform_into_standard_form()
Definition: polynomial.cc:112
void scale(const number factor)
Definition: polynomial.cc:166
Polynomial< number > & operator-=(const Polynomial< number > &p)
Definition: polynomial.cc:311
std::vector< number > lagrange_support_points
Definition: polynomial.h:313
void shift(const number2 offset)
Definition: polynomial.cc:440
void print(std::ostream &out) const
Definition: polynomial.cc:515
void serialize(Archive &ar, const unsigned int version)
Definition: polynomial.h:1023
static void multiply(std::vector< number > &coefficients, const number factor)
Definition: polynomial.cc:191
void value(const Number2 x, const unsigned int n_derivatives, Number2 *values) const
Definition: polynomial.h:846
Polynomial< number > & operator*=(const double s)
Definition: polynomial.cc:204
virtual std::size_t memory_consumption() const
Definition: polynomial.cc:534
unsigned int degree() const
Definition: polynomial.h:799
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 2 > second
Definition: grid_out.cc:4604
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcEmptyObject()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcMessage(std::string arg1)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
Number jacobi_polynomial_value(const unsigned int degree, const int alpha, const int beta, const Number x)
Definition: polynomial.h:1037
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 > > &points)
Definition: polynomial.cc:702
std::vector< Number > jacobi_polynomial_roots(const unsigned int degree, const int alpha, const int beta)
Definition: polynomial.h:1078
static constexpr double PI
Definition: numbers.h:233
static const unsigned int invalid_unsigned_int
Definition: types.h:201
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)