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