Reference documentation for deal.II version GIT 574c7c8486 2023-09-22 10:20: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 <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] = weight;
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 derivatives by k! to transform the product
901  // p_n = p^(n)(x)/k! into the actual form of the derivative
902  {
903  number k_factorial = 2;
904  for (unsigned int k = 2; k <= n_derivatives; ++k)
905  {
906  for (unsigned int e = 0; e < n_entries; ++e)
907  values[k][e] *= k_factorial;
908  k_factorial *= static_cast<number>(k + 1);
909  }
910  }
911  break;
912 
913  // manually implement case 0 (values only), case 1 (value + first
914  // derivative), and case 2 (up to second derivative) since they
915  // might be called often. then, we can unroll the inner loop and
916  // keep the temporary results as local variables to help the
917  // compiler with the pointer aliasing analysis.
918  case 0:
919  {
920  std::array<Number2, n_entries> value;
921  for (unsigned int e = 0; e < n_entries; ++e)
922  value[e] = weight;
923  for (unsigned int i = 0; i < n_supp; ++i)
924  for (unsigned int e = 0; e < n_entries; ++e)
925  value[e] *= (x[e] - lagrange_support_points[i]);
926 
927  for (unsigned int e = 0; e < n_entries; ++e)
928  values[0][e] = value[e];
929  break;
930  }
931 
932  case 1:
933  {
934  std::array<Number2, n_entries> value, derivative = {};
935  for (unsigned int e = 0; e < n_entries; ++e)
936  value[e] = weight;
937  for (unsigned int i = 0; i < n_supp; ++i)
938  for (unsigned int e = 0; e < n_entries; ++e)
939  {
940  const Number2 v = x[e] - lagrange_support_points[i];
941  derivative[e] = derivative[e] * v + value[e];
942  value[e] *= v;
943  }
944 
945  for (unsigned int e = 0; e < n_entries; ++e)
946  {
947  values[0][e] = value[e];
948  values[1][e] = derivative[e];
949  }
950  break;
951  }
952 
953  case 2:
954  {
955  std::array<Number2, n_entries> value, derivative = {},
956  second = {};
957  for (unsigned int e = 0; e < n_entries; ++e)
958  value[e] = weight;
959  for (unsigned int i = 0; i < n_supp; ++i)
960  for (unsigned int e = 0; e < n_entries; ++e)
961  {
962  const Number2 v = x[e] - lagrange_support_points[i];
963  second[e] = second[e] * v + derivative[e];
964  derivative[e] = derivative[e] * v + value[e];
965  value[e] *= v;
966  }
967 
968  for (unsigned int e = 0; e < n_entries; ++e)
969  {
970  values[0][e] = value[e];
971  values[1][e] = derivative[e];
972  values[2][e] = static_cast<number>(2) * second[e];
973  }
974  break;
975  }
976  }
977  return;
978  }
979 
980  Assert(coefficients.size() > 0, ExcEmptyObject());
981 
982  // if derivatives are needed, then do it properly by the full
983  // Horner scheme
984  const unsigned int m = coefficients.size();
985  std::vector<std::array<Number2, n_entries>> a(coefficients.size());
986  for (unsigned int i = 0; i < coefficients.size(); ++i)
987  for (unsigned int e = 0; e < n_entries; ++e)
988  a[i][e] = coefficients[i];
989 
990  unsigned int j_factorial = 1;
991 
992  // loop over all requested derivatives. note that derivatives @p{j>m} are
993  // necessarily zero, as they differentiate the polynomial more often than
994  // the highest power is
995  const unsigned int min_valuessize_m = std::min(n_derivatives + 1, m);
996  for (unsigned int j = 0; j < min_valuessize_m; ++j)
997  {
998  for (int k = m - 2; k >= static_cast<int>(j); --k)
999  for (unsigned int e = 0; e < n_entries; ++e)
1000  a[k][e] += x[e] * a[k + 1][e];
1001  for (unsigned int e = 0; e < n_entries; ++e)
1002  values[j][e] = static_cast<number>(j_factorial) * a[j][e];
1003 
1004  j_factorial *= j + 1;
1005  }
1006 
1007  // fill higher derivatives by zero
1008  for (unsigned int j = min_valuessize_m; j <= n_derivatives; ++j)
1009  for (unsigned int e = 0; e < n_entries; ++e)
1010  values[j][e] = 0.;
1011  }
1012 
1013 
1014 
1015  template <typename number>
1016  template <class Archive>
1017  inline void
1018  Polynomial<number>::serialize(Archive &ar, const unsigned int)
1019  {
1020  // forward to serialization function in the base class.
1021  ar &static_cast<Subscriptor &>(*this);
1022  ar &coefficients;
1023  ar &in_lagrange_product_form;
1024  ar &lagrange_support_points;
1025  ar &lagrange_weight;
1026  }
1027 
1028 
1029 
1030  template <typename Number>
1031  Number
1032  jacobi_polynomial_value(const unsigned int degree,
1033  const int alpha,
1034  const int beta,
1035  const Number x)
1036  {
1037  Assert(alpha >= 0 && beta >= 0,
1038  ExcNotImplemented("Negative alpha/beta coefficients not supported"));
1039  // the Jacobi polynomial is evaluated using a recursion formula.
1040  Number p0, p1;
1041 
1042  // The recursion formula is defined for the interval [-1, 1], so rescale
1043  // to that interval here
1044  const Number xeval = Number(-1) + 2. * x;
1045 
1046  // initial values P_0(x), P_1(x):
1047  p0 = 1.0;
1048  if (degree == 0)
1049  return p0;
1050  p1 = ((alpha + beta + 2) * xeval + (alpha - beta)) / 2;
1051  if (degree == 1)
1052  return p1;
1053 
1054  for (unsigned int i = 1; i < degree; ++i)
1055  {
1056  const Number v = 2 * i + (alpha + beta);
1057  const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
1058  const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
1059  const Number a3 = v * (v + 1) * (v + 2);
1060  const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
1061 
1062  const Number pn = ((a2 + a3 * xeval) * p1 - a4 * p0) / a1;
1063  p0 = p1;
1064  p1 = pn;
1065  }
1066  return p1;
1067  }
1068 
1069 
1070 
1071  template <typename Number>
1072  std::vector<Number>
1073  jacobi_polynomial_roots(const unsigned int degree,
1074  const int alpha,
1075  const int beta)
1076  {
1077  std::vector<Number> x(degree, 0.5);
1078 
1079  // compute zeros with a Newton algorithm.
1080 
1081  // Set tolerance. For long double we might not always get the additional
1082  // precision in a run time environment (e.g. with valgrind), so we must
1083  // limit the tolerance to double. Since we do a Newton iteration, doing
1084  // one more iteration after the residual has indicated convergence will be
1085  // enough for all number types due to the quadratic convergence of
1086  // Newton's method
1087 
1088  const Number tolerance =
1089  4 * std::max(static_cast<Number>(std::numeric_limits<double>::epsilon()),
1091 
1092  // The following implementation follows closely the one given in the
1093  // appendix of the book by Karniadakis and Sherwin: Spectral/hp element
1094  // methods for computational fluid dynamics (Oxford University Press,
1095  // 2005)
1096 
1097  // If symmetric, we only need to compute the half of points
1098  const unsigned int n_points = (alpha == beta ? degree / 2 : degree);
1099  for (unsigned int k = 0; k < n_points; ++k)
1100  {
1101  // we take the zeros of the Chebyshev polynomial (alpha=beta=-0.5) as
1102  // initial values, corrected by the initial value
1103  Number r = 0.5 - 0.5 * std::cos(static_cast<Number>(2 * k + 1) /
1104  (2 * degree) * numbers::PI);
1105  if (k > 0)
1106  r = (r + x[k - 1]) / 2;
1107 
1108  unsigned int converged = numbers::invalid_unsigned_int;
1109  for (unsigned int it = 1; it < 1000; ++it)
1110  {
1111  Number s = 0.;
1112  for (unsigned int i = 0; i < k; ++i)
1113  s += 1. / (r - x[i]);
1114 
1115  // derivative of P_n^{alpha,beta}, rescaled to [0, 1]
1116  const Number J_x =
1117  (alpha + beta + degree + 1) *
1118  jacobi_polynomial_value(degree - 1, alpha + 1, beta + 1, r);
1119 
1120  // value of P_n^{alpha,beta}
1121  const Number f = jacobi_polynomial_value(degree, alpha, beta, r);
1122  const Number delta = f / (f * s - J_x);
1123  r += delta;
1124  if (converged == numbers::invalid_unsigned_int &&
1125  std::abs(delta) < tolerance)
1126  converged = it;
1127 
1128  // do one more iteration to ensure accuracy also for tighter
1129  // types than double (e.g. long double)
1130  if (it == converged + 1)
1131  break;
1132  }
1133 
1135  ExcMessage("Newton iteration for zero of Jacobi polynomial "
1136  "did not converge."));
1137 
1138  x[k] = r;
1139  }
1140 
1141  // in case we assumed symmetry, fill up the missing values
1142  for (unsigned int k = n_points; k < degree; ++k)
1143  x[k] = 1.0 - x[degree - k - 1];
1144 
1145  return x;
1146  }
1147 
1148 } // namespace Polynomials
1150 
1151 #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:1018
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: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:1616
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:1032
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:1073
static constexpr double PI
Definition: numbers.h:260
static const unsigned int invalid_unsigned_int
Definition: types.h:213