Reference documentation for deal.II version GIT 9042b9283b 2023-12-02 14: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 - 2023 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 <limits>
29 #include <memory>
30 #include <shared_mutex>
31 #include <vector>
32 
34 
44 namespace Polynomials
45 {
65  template <typename number>
66  class Polynomial : public Subscriptor
67  {
68  public:
77  Polynomial(const std::vector<number> &coefficients);
78 
82  Polynomial(const unsigned int n);
83 
91  Polynomial(const std::vector<Point<1>> &lagrange_support_points,
92  const unsigned int evaluation_point);
93 
98 
108  number
109  value(const number x) const;
110 
121  void
122  value(const number x, std::vector<number> &values) const;
123 
142  template <typename Number2>
143  void
144  value(const Number2 x,
145  const unsigned int n_derivatives,
146  Number2 *values) const;
147 
160  template <std::size_t n_entries, typename Number2>
161  void
162  values_of_array(const std::array<Number2, n_entries> &points,
163  const unsigned int n_derivatives,
164  std::array<Number2, n_entries> *values) const;
165 
171  unsigned int
172  degree() const;
173 
181  void
182  scale(const number factor);
183 
199  template <typename number2>
200  void
201  shift(const number2 offset);
202 
207  derivative() const;
208 
214  primitive() const;
215 
220  operator*=(const double s);
221 
226  operator*=(const Polynomial<number> &p);
227 
232  operator+=(const Polynomial<number> &p);
233 
238  operator-=(const Polynomial<number> &p);
239 
243  bool
244  operator==(const Polynomial<number> &p) const;
245 
249  void
250  print(std::ostream &out) const;
251 
257  template <class Archive>
258  void
259  serialize(Archive &ar, const unsigned int version);
260 
264  virtual std::size_t
265  memory_consumption() const;
266 
267  protected:
271  static void
272  scale(std::vector<number> &coefficients, const number factor);
273 
277  template <typename number2>
278  static void
279  shift(std::vector<number> &coefficients, const number2 shift);
280 
284  static void
285  multiply(std::vector<number> &coefficients, const number factor);
286 
292  void
294 
303  std::vector<number> coefficients;
304 
310 
315  std::vector<number> lagrange_support_points;
316 
322  };
323 
324 
329  template <typename number>
330  class Monomial : public Polynomial<number>
331  {
332  public:
337  Monomial(const unsigned int n, const double coefficient = 1.);
338 
345  static std::vector<Polynomial<number>>
346  generate_complete_basis(const unsigned int degree);
347 
348  private:
352  static std::vector<number>
353  make_vector(unsigned int n, const double coefficient);
354  };
355 
356 
372  class LagrangeEquidistant : public Polynomial<double>
373  {
374  public:
380  LagrangeEquidistant(const unsigned int n, const unsigned int support_point);
381 
390  static std::vector<Polynomial<double>>
391  generate_complete_basis(const unsigned int degree);
392 
393  private:
398  static void
399  compute_coefficients(const unsigned int n,
400  const unsigned int support_point,
401  std::vector<double> &a);
402  };
403 
404 
405 
412  std::vector<Polynomial<double>>
413  generate_complete_Lagrange_basis(const std::vector<Point<1>> &points);
414 
415 
416 
430  class Legendre : public Polynomial<double>
431  {
432  public:
436  Legendre(const unsigned int p);
437 
444  static std::vector<Polynomial<double>>
445  generate_complete_basis(const unsigned int degree);
446  };
447 
467  class Lobatto : public Polynomial<double>
468  {
469  public:
474  Lobatto(const unsigned int p = 0);
475 
480  static std::vector<Polynomial<double>>
481  generate_complete_basis(const unsigned int p);
482 
483  private:
487  std::vector<double>
488  compute_coefficients(const unsigned int p);
489  };
490 
491 
492 
530  class Hierarchical : public Polynomial<double>
531  {
532  public:
537  Hierarchical(const unsigned int p);
538 
549  static std::vector<Polynomial<double>>
550  generate_complete_basis(const unsigned int degree);
551 
552  private:
556  static void
557  compute_coefficients(const unsigned int p);
558 
563  static const std::vector<double> &
564  get_coefficients(const unsigned int p);
565 
574  static std::vector<std::unique_ptr<const std::vector<double>>>
576 
581  static std::shared_mutex coefficients_lock;
582  };
583 
584 
585 
613  class HermiteInterpolation : public Polynomial<double>
614  {
615  public:
620  HermiteInterpolation(const unsigned int p);
621 
627  static std::vector<Polynomial<double>>
628  generate_complete_basis(const unsigned int p);
629  };
630 
631 
632 
734  class HermiteLikeInterpolation : public Polynomial<double>
735  {
736  public:
741  HermiteLikeInterpolation(const unsigned int degree,
742  const unsigned int index);
743 
748  static std::vector<Polynomial<double>>
749  generate_complete_basis(const unsigned int degree);
750  };
751 
752 
753 
754  /*
755  * Evaluate a Jacobi polynomial @f$ P_n^{\alpha, \beta}(x) @f$ specified by the
756  * parameters @p alpha, @p beta, @p n, where @p n is the degree of the
757  * Jacobi polynomial.
758  *
759  * @note The Jacobi polynomials are not orthonormal and are defined on the
760  * unit interval @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often
761  * used in literature. @p x is the point of evaluation.
762  */
763  template <typename Number>
764  Number
765  jacobi_polynomial_value(const unsigned int degree,
766  const int alpha,
767  const int beta,
768  const Number x);
769 
770 
783  template <typename Number>
784  std::vector<Number>
785  jacobi_polynomial_roots(const unsigned int degree,
786  const int alpha,
787  const int beta);
788 } // namespace Polynomials
789 
790 
793 /* -------------------------- inline functions --------------------- */
794 
795 namespace Polynomials
796 {
797  template <typename number>
799  : in_lagrange_product_form(false)
800  , lagrange_weight(1.)
801  {}
802 
803 
804 
805  template <typename number>
806  inline unsigned int
808  {
809  if (in_lagrange_product_form == true)
810  {
811  return lagrange_support_points.size();
812  }
813  else
814  {
815  Assert(coefficients.size() > 0, ExcEmptyObject());
816  return coefficients.size() - 1;
817  }
818  }
819 
820 
821 
822  template <typename number>
823  inline number
824  Polynomial<number>::value(const number x) const
825  {
826  if (in_lagrange_product_form == false)
827  {
828  Assert(coefficients.size() > 0, ExcEmptyObject());
829 
830  // Horner scheme
831  const unsigned int m = coefficients.size();
832  number value = coefficients.back();
833  for (int k = m - 2; k >= 0; --k)
834  value = value * x + coefficients[k];
835  return value;
836  }
837  else
838  {
839  // direct evaluation of Lagrange polynomial
840  const unsigned int m = lagrange_support_points.size();
841  number value = 1.;
842  for (unsigned int j = 0; j < m; ++j)
843  value *= x - lagrange_support_points[j];
844  value *= lagrange_weight;
845  return value;
846  }
847  }
848 
849 
850 
851  template <typename number>
852  template <typename Number2>
853  inline void
854  Polynomial<number>::value(const Number2 x,
855  const unsigned int n_derivatives,
856  Number2 *values) const
857  {
858  values_of_array(std::array<Number2, 1ul>{{x}},
859  n_derivatives,
860  reinterpret_cast<std::array<Number2, 1ul> *>(values));
861  }
862 
863 
864 
865  template <typename number>
866  template <std::size_t n_entries, typename Number2>
867  inline void
869  const std::array<Number2, n_entries> &x,
870  const unsigned int n_derivatives,
871  std::array<Number2, n_entries> *values) const
872  {
873  // evaluate Lagrange polynomial and derivatives
874  if (in_lagrange_product_form == true)
875  {
876  // to compute the value and all derivatives of a polynomial of the
877  // form (x-x_1)*(x-x_2)*...*(x-x_n), expand the derivatives like
878  // automatic differentiation does.
879  const unsigned int n_supp = lagrange_support_points.size();
880  const number weight = lagrange_weight;
881  switch (n_derivatives)
882  {
883  default:
884  for (unsigned int e = 0; e < n_entries; ++e)
885  values[0][e] = weight;
886  for (unsigned int k = 1; k <= n_derivatives; ++k)
887  for (unsigned int e = 0; e < n_entries; ++e)
888  values[k][e] = 0.;
889  for (unsigned int i = 0; i < n_supp; ++i)
890  {
891  std::array<Number2, n_entries> v = x;
892  for (unsigned int e = 0; e < n_entries; ++e)
893  v[e] -= lagrange_support_points[i];
894 
895  // multiply by (x-x_i) and compute action on all derivatives,
896  // too (inspired from automatic differentiation: implement the
897  // product rule for the old value and the new variable 'v',
898  // i.e., expand value v and derivative one). since we reuse a
899  // value from the next lower derivative from the steps before,
900  // need to start from the highest derivative
901  for (unsigned int k = n_derivatives; k > 0; --k)
902  for (unsigned int e = 0; e < n_entries; ++e)
903  values[k][e] = (values[k][e] * v[e] + values[k - 1][e]);
904  for (unsigned int e = 0; e < n_entries; ++e)
905  values[0][e] *= v[e];
906  }
907  // finally, multiply derivatives by k! to transform the product
908  // p_n = p^(n)(x)/k! into the actual form of the derivative
909  {
910  number k_factorial = 2;
911  for (unsigned int k = 2; k <= n_derivatives; ++k)
912  {
913  for (unsigned int e = 0; e < n_entries; ++e)
914  values[k][e] *= k_factorial;
915  k_factorial *= static_cast<number>(k + 1);
916  }
917  }
918  break;
919 
920  // manually implement case 0 (values only), case 1 (value + first
921  // derivative), and case 2 (up to second derivative) since they
922  // might be called often. then, we can unroll the inner loop and
923  // keep the temporary results as local variables to help the
924  // compiler with the pointer aliasing analysis.
925  case 0:
926  {
927  std::array<Number2, n_entries> value;
928  for (unsigned int e = 0; e < n_entries; ++e)
929  value[e] = weight;
930  for (unsigned int i = 0; i < n_supp; ++i)
931  for (unsigned int e = 0; e < n_entries; ++e)
932  value[e] *= (x[e] - lagrange_support_points[i]);
933 
934  for (unsigned int e = 0; e < n_entries; ++e)
935  values[0][e] = value[e];
936  break;
937  }
938 
939  case 1:
940  {
941  std::array<Number2, n_entries> value, derivative = {};
942  for (unsigned int e = 0; e < n_entries; ++e)
943  value[e] = weight;
944  for (unsigned int i = 0; i < n_supp; ++i)
945  for (unsigned int e = 0; e < n_entries; ++e)
946  {
947  const Number2 v = x[e] - lagrange_support_points[i];
948  derivative[e] = derivative[e] * v + value[e];
949  value[e] *= v;
950  }
951 
952  for (unsigned int e = 0; e < n_entries; ++e)
953  {
954  values[0][e] = value[e];
955  values[1][e] = derivative[e];
956  }
957  break;
958  }
959 
960  case 2:
961  {
962  std::array<Number2, n_entries> value, derivative = {},
963  second = {};
964  for (unsigned int e = 0; e < n_entries; ++e)
965  value[e] = weight;
966  for (unsigned int i = 0; i < n_supp; ++i)
967  for (unsigned int e = 0; e < n_entries; ++e)
968  {
969  const Number2 v = x[e] - lagrange_support_points[i];
970  second[e] = second[e] * v + derivative[e];
971  derivative[e] = derivative[e] * v + value[e];
972  value[e] *= v;
973  }
974 
975  for (unsigned int e = 0; e < n_entries; ++e)
976  {
977  values[0][e] = value[e];
978  values[1][e] = derivative[e];
979  values[2][e] = static_cast<number>(2) * second[e];
980  }
981  break;
982  }
983  }
984  return;
985  }
986 
987  Assert(coefficients.size() > 0, ExcEmptyObject());
988 
989  // if derivatives are needed, then do it properly by the full
990  // Horner scheme
991  const unsigned int m = coefficients.size();
992  std::vector<std::array<Number2, n_entries>> a(coefficients.size());
993  for (unsigned int i = 0; i < coefficients.size(); ++i)
994  for (unsigned int e = 0; e < n_entries; ++e)
995  a[i][e] = coefficients[i];
996 
997  unsigned int j_factorial = 1;
998 
999  // loop over all requested derivatives. note that derivatives @p{j>m} are
1000  // necessarily zero, as they differentiate the polynomial more often than
1001  // the highest power is
1002  const unsigned int min_valuessize_m = std::min(n_derivatives + 1, m);
1003  for (unsigned int j = 0; j < min_valuessize_m; ++j)
1004  {
1005  for (int k = m - 2; k >= static_cast<int>(j); --k)
1006  for (unsigned int e = 0; e < n_entries; ++e)
1007  a[k][e] += x[e] * a[k + 1][e];
1008  for (unsigned int e = 0; e < n_entries; ++e)
1009  values[j][e] = static_cast<number>(j_factorial) * a[j][e];
1010 
1011  j_factorial *= j + 1;
1012  }
1013 
1014  // fill higher derivatives by zero
1015  for (unsigned int j = min_valuessize_m; j <= n_derivatives; ++j)
1016  for (unsigned int e = 0; e < n_entries; ++e)
1017  values[j][e] = 0.;
1018  }
1019 
1020 
1021 
1022  template <typename number>
1023  template <class Archive>
1024  inline void
1025  Polynomial<number>::serialize(Archive &ar, const unsigned int)
1026  {
1027  // forward to serialization function in the base class.
1028  ar &static_cast<Subscriptor &>(*this);
1029  ar &coefficients;
1030  ar &in_lagrange_product_form;
1031  ar &lagrange_support_points;
1032  ar &lagrange_weight;
1033  }
1034 
1035 
1036 
1037  template <typename Number>
1038  Number
1039  jacobi_polynomial_value(const unsigned int degree,
1040  const int alpha,
1041  const int beta,
1042  const Number x)
1043  {
1044  Assert(alpha >= 0 && beta >= 0,
1045  ExcNotImplemented("Negative alpha/beta coefficients not supported"));
1046  // the Jacobi polynomial is evaluated using a recursion formula.
1047  Number p0, p1;
1048 
1049  // The recursion formula is defined for the interval [-1, 1], so rescale
1050  // to that interval here
1051  const Number xeval = Number(-1) + 2. * x;
1052 
1053  // initial values P_0(x), P_1(x):
1054  p0 = 1.0;
1055  if (degree == 0)
1056  return p0;
1057  p1 = ((alpha + beta + 2) * xeval + (alpha - beta)) / 2;
1058  if (degree == 1)
1059  return p1;
1060 
1061  for (unsigned int i = 1; i < degree; ++i)
1062  {
1063  const Number v = 2 * i + (alpha + beta);
1064  const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
1065  const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
1066  const Number a3 = v * (v + 1) * (v + 2);
1067  const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
1068 
1069  const Number pn = ((a2 + a3 * xeval) * p1 - a4 * p0) / a1;
1070  p0 = p1;
1071  p1 = pn;
1072  }
1073  return p1;
1074  }
1075 
1076 
1077 
1078  template <typename Number>
1079  std::vector<Number>
1080  jacobi_polynomial_roots(const unsigned int degree,
1081  const int alpha,
1082  const int beta)
1083  {
1084  std::vector<Number> x(degree, 0.5);
1085 
1086  // compute zeros with a Newton algorithm.
1087 
1088  // Set tolerance. For long double we might not always get the additional
1089  // precision in a run time environment (e.g. with valgrind), so we must
1090  // limit the tolerance to double. Since we do a Newton iteration, doing
1091  // one more iteration after the residual has indicated convergence will be
1092  // enough for all number types due to the quadratic convergence of
1093  // Newton's method
1094 
1095  const Number tolerance =
1096  4 * std::max(static_cast<Number>(std::numeric_limits<double>::epsilon()),
1098 
1099  // The following implementation follows closely the one given in the
1100  // appendix of the book by Karniadakis and Sherwin: Spectral/hp element
1101  // methods for computational fluid dynamics (Oxford University Press,
1102  // 2005)
1103 
1104  // If symmetric, we only need to compute the half of points
1105  const unsigned int n_points = (alpha == beta ? degree / 2 : degree);
1106  for (unsigned int k = 0; k < n_points; ++k)
1107  {
1108  // we take the zeros of the Chebyshev polynomial (alpha=beta=-0.5) as
1109  // initial values, corrected by the initial value
1110  Number r = 0.5 - 0.5 * std::cos(static_cast<Number>(2 * k + 1) /
1111  (2 * degree) * numbers::PI);
1112  if (k > 0)
1113  r = (r + x[k - 1]) / 2;
1114 
1115  unsigned int converged = numbers::invalid_unsigned_int;
1116  for (unsigned int it = 1; it < 1000; ++it)
1117  {
1118  Number s = 0.;
1119  for (unsigned int i = 0; i < k; ++i)
1120  s += 1. / (r - x[i]);
1121 
1122  // derivative of P_n^{alpha,beta}, rescaled to [0, 1]
1123  const Number J_x =
1124  (alpha + beta + degree + 1) *
1125  jacobi_polynomial_value(degree - 1, alpha + 1, beta + 1, r);
1126 
1127  // value of P_n^{alpha,beta}
1128  const Number f = jacobi_polynomial_value(degree, alpha, beta, r);
1129  const Number delta = f / (f * s - J_x);
1130  r += delta;
1131  if (converged == numbers::invalid_unsigned_int &&
1132  std::abs(delta) < tolerance)
1133  converged = it;
1134 
1135  // do one more iteration to ensure accuracy also for tighter
1136  // types than double (e.g. long double)
1137  if (it == converged + 1)
1138  break;
1139  }
1140 
1142  ExcMessage("Newton iteration for zero of Jacobi polynomial "
1143  "did not converge."));
1144 
1145  x[k] = r;
1146  }
1147 
1148  // in case we assumed symmetry, fill up the missing values
1149  for (unsigned int k = n_points; k < degree; ++k)
1150  x[k] = 1.0 - x[degree - k - 1];
1151 
1152  return x;
1153  }
1154 
1155 } // namespace Polynomials
1157 
1158 #endif
HermiteInterpolation(const unsigned int p)
Definition: polynomial.cc:1049
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
Definition: polynomial.cc:1101
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:1397
HermiteLikeInterpolation(const unsigned int degree, const unsigned int index)
Definition: polynomial.cc:1167
static std::vector< std::unique_ptr< const std::vector< double > > > recursive_coefficients
Definition: polynomial.h:575
Hierarchical(const unsigned int p)
Definition: polynomial.cc:857
static void compute_coefficients(const unsigned int p)
Definition: polynomial.cc:864
static const std::vector< double > & get_coefficients(const unsigned int p)
Definition: polynomial.cc:1006
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:1021
static std::shared_mutex coefficients_lock
Definition: polynomial.h:581
static void compute_coefficients(const unsigned int n, const unsigned int support_point, std::vector< double > &a)
Definition: polynomial.cc:606
LagrangeEquidistant(const unsigned int n, const unsigned int support_point)
Definition: polynomial.cc:582
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:664
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:730
Legendre(const unsigned int p)
Definition: polynomial.cc:703
std::vector< double > compute_coefficients(const unsigned int p)
Definition: polynomial.cc:749
Lobatto(const unsigned int p=0)
Definition: polynomial.cc:744
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
Definition: polynomial.cc:835
static std::vector< Polynomial< number > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:551
Monomial(const unsigned int n, const double coefficient=1.)
Definition: polynomial.cc:543
static std::vector< number > make_vector(unsigned int n, const double coefficient)
Definition: polynomial.cc:533
number value(const number x) const
Definition: polynomial.h:824
bool operator==(const Polynomial< number > &p) const
Definition: polynomial.cc:332
std::vector< number > coefficients
Definition: polynomial.h:303
Polynomial< number > primitive() const
Definition: polynomial.cc:473
Polynomial< number > & operator+=(const Polynomial< number > &p)
Definition: polynomial.cc:254
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:868
Polynomial< number > derivative() const
Definition: polynomial.cc:444
void transform_into_standard_form()
Definition: polynomial.cc:97
void scale(const number factor)
Definition: polynomial.cc:151
Polynomial< number > & operator-=(const Polynomial< number > &p)
Definition: polynomial.cc:296
std::vector< number > lagrange_support_points
Definition: polynomial.h:315
void shift(const number2 offset)
Definition: polynomial.cc:425
void print(std::ostream &out) const
Definition: polynomial.cc:500
void serialize(Archive &ar, const unsigned int version)
Definition: polynomial.h:1025
static void multiply(std::vector< number > &coefficients, const number factor)
Definition: polynomial.cc:176
void value(const Number2 x, const unsigned int n_derivatives, Number2 *values) const
Definition: polynomial.h:854
Polynomial< number > & operator*=(const double s)
Definition: polynomial.cc:189
virtual std::size_t memory_consumption() const
Definition: polynomial.cc:519
unsigned int degree() const
Definition: polynomial.h:807
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 2 > second
Definition: grid_out.cc:4615
#define Assert(cond, exc)
Definition: exceptions.h:1631
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:1039
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
Definition: polynomial.cc:687
std::vector< Number > jacobi_polynomial_roots(const unsigned int degree, const int alpha, const int beta)
Definition: polynomial.h:1080
static constexpr double PI
Definition: numbers.h:260
static const unsigned int invalid_unsigned_int
Definition: types.h:221