Reference documentation for deal.II version GIT d7aca55de5 2022-08-10 12:50:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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
9 // Public License as published by the Free Software Foundation; either
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 
42 namespace 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 
89  Polynomial(const std::vector<Point<1>> &lagrange_support_points,
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 
224  operator*=(const Polynomial<number> &p);
225 
230  operator+=(const Polynomial<number> &p);
231 
236  operator-=(const Polynomial<number> &p);
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 
726  class HermiteLikeInterpolation : public Polynomial<double>
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 
787 namespace 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
846  Polynomial<number>::value(const Number2 x,
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:4605
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcEmptyObject()
static ::ExceptionBase & ExcMessage(std::string arg1)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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:248
static const unsigned int invalid_unsigned_int
Definition: types.h:201