Reference documentation for deal.II version Git 4fbb5374f0 2021-01-22 15:02:13 -0500
\(\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 - 2019 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 <memory>
28 #include <vector>
29 
31 
41 namespace Polynomials
42 {
59  template <typename number>
60  class Polynomial : public Subscriptor
61  {
62  public:
71  Polynomial(const std::vector<number> &coefficients);
72 
76  Polynomial(const unsigned int n);
77 
85  Polynomial(const std::vector<Point<1>> &lagrange_support_points,
86  const unsigned int evaluation_point);
87 
91  Polynomial();
92 
102  number
103  value(const number x) const;
104 
115  void
116  value(const number x, std::vector<number> &values) const;
117 
136  template <typename Number2>
137  void
138  value(const Number2 x,
139  const unsigned int n_derivatives,
140  Number2 * values) const;
141 
147  unsigned int
148  degree() const;
149 
157  void
158  scale(const number factor);
159 
175  template <typename number2>
176  void
177  shift(const number2 offset);
178 
183  derivative() const;
184 
190  primitive() const;
191 
196  operator*=(const double s);
197 
202  operator*=(const Polynomial<number> &p);
203 
208  operator+=(const Polynomial<number> &p);
209 
214  operator-=(const Polynomial<number> &p);
215 
219  bool
220  operator==(const Polynomial<number> &p) const;
221 
225  void
226  print(std::ostream &out) const;
227 
233  template <class Archive>
234  void
235  serialize(Archive &ar, const unsigned int version);
236 
240  virtual std::size_t
241  memory_consumption() const;
242 
243  protected:
247  static void
248  scale(std::vector<number> &coefficients, const number factor);
249 
253  template <typename number2>
254  static void
255  shift(std::vector<number> &coefficients, const number2 shift);
256 
260  static void
261  multiply(std::vector<number> &coefficients, const number factor);
262 
268  void
270 
279  std::vector<number> coefficients;
280 
286 
291  std::vector<number> lagrange_support_points;
292 
298  };
299 
300 
305  template <typename number>
306  class Monomial : public Polynomial<number>
307  {
308  public:
313  Monomial(const unsigned int n, const double coefficient = 1.);
314 
321  static std::vector<Polynomial<number>>
322  generate_complete_basis(const unsigned int degree);
323 
324  private:
328  static std::vector<number>
329  make_vector(unsigned int n, const double coefficient);
330  };
331 
332 
348  class LagrangeEquidistant : public Polynomial<double>
349  {
350  public:
356  LagrangeEquidistant(const unsigned int n, const unsigned int support_point);
357 
366  static std::vector<Polynomial<double>>
367  generate_complete_basis(const unsigned int degree);
368 
369  private:
374  static void
375  compute_coefficients(const unsigned int n,
376  const unsigned int support_point,
377  std::vector<double> &a);
378  };
379 
380 
381 
388  std::vector<Polynomial<double>>
389  generate_complete_Lagrange_basis(const std::vector<Point<1>> &points);
390 
391 
392 
406  class Legendre : public Polynomial<double>
407  {
408  public:
412  Legendre(const unsigned int p);
413 
420  static std::vector<Polynomial<double>>
421  generate_complete_basis(const unsigned int degree);
422  };
423 
443  class Lobatto : public Polynomial<double>
444  {
445  public:
450  Lobatto(const unsigned int p = 0);
451 
456  static std::vector<Polynomial<double>>
457  generate_complete_basis(const unsigned int p);
458 
459  private:
463  std::vector<double>
464  compute_coefficients(const unsigned int p);
465  };
466 
467 
468 
506  class Hierarchical : public Polynomial<double>
507  {
508  public:
513  Hierarchical(const unsigned int p);
514 
525  static std::vector<Polynomial<double>>
526  generate_complete_basis(const unsigned int degree);
527 
528  private:
532  static void
533  compute_coefficients(const unsigned int p);
534 
539  static const std::vector<double> &
540  get_coefficients(const unsigned int p);
541 
550  static std::vector<std::unique_ptr<const std::vector<double>>>
552  };
553 
554 
555 
583  class HermiteInterpolation : public Polynomial<double>
584  {
585  public:
590  HermiteInterpolation(const unsigned int p);
591 
597  static std::vector<Polynomial<double>>
598  generate_complete_basis(const unsigned int p);
599  };
600 
601 
602 
704  class HermiteLikeInterpolation : public Polynomial<double>
705  {
706  public:
711  HermiteLikeInterpolation(const unsigned int degree,
712  const unsigned int index);
713 
718  static std::vector<Polynomial<double>>
719  generate_complete_basis(const unsigned int degree);
720  };
721 
722 
723 
724  /*
725  * Evaluate a Jacobi polynomial @f$ P_n^{\alpha, \beta}(x) @f$ specified by the
726  * parameters @p alpha, @p beta, @p n, where @p n is the degree of the
727  * Jacobi polynomial.
728  *
729  * @note The Jacobi polynomials are not orthonormal and are defined on the
730  * unit interval @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often
731  * used in literature. @p x is the point of evaluation.
732  */
733  template <typename Number>
734  Number
735  jacobi_polynomial_value(const unsigned int degree,
736  const int alpha,
737  const int beta,
738  const Number x);
739 
740 
753  template <typename Number>
754  std::vector<Number>
755  jacobi_polynomial_roots(const unsigned int degree,
756  const int alpha,
757  const int beta);
758 } // namespace Polynomials
759 
760 
763 /* -------------------------- inline functions --------------------- */
764 
765 namespace Polynomials
766 {
767  template <typename number>
769  : in_lagrange_product_form(false)
770  , lagrange_weight(1.)
771  {}
772 
773 
774 
775  template <typename number>
776  inline unsigned int
778  {
779  if (in_lagrange_product_form == true)
780  {
781  return lagrange_support_points.size();
782  }
783  else
784  {
785  Assert(coefficients.size() > 0, ExcEmptyObject());
786  return coefficients.size() - 1;
787  }
788  }
789 
790 
791 
792  template <typename number>
793  inline number
794  Polynomial<number>::value(const number x) const
795  {
796  if (in_lagrange_product_form == false)
797  {
798  Assert(coefficients.size() > 0, ExcEmptyObject());
799 
800  // Horner scheme
801  const unsigned int m = coefficients.size();
802  number value = coefficients.back();
803  for (int k = m - 2; k >= 0; --k)
804  value = value * x + coefficients[k];
805  return value;
806  }
807  else
808  {
809  // direct evaluation of Lagrange polynomial
810  const unsigned int m = lagrange_support_points.size();
811  number value = 1.;
812  for (unsigned int j = 0; j < m; ++j)
813  value *= x - lagrange_support_points[j];
814  value *= lagrange_weight;
815  return value;
816  }
817  }
818 
819 
820 
821  template <typename number>
822  template <typename Number2>
823  inline void
824  Polynomial<number>::value(const Number2 x,
825  const unsigned int n_derivatives,
826  Number2 * values) const
827  {
828  // evaluate Lagrange polynomial and derivatives
829  if (in_lagrange_product_form == true)
830  {
831  // to compute the value and all derivatives of a polynomial of the
832  // form (x-x_1)*(x-x_2)*...*(x-x_n), expand the derivatives like
833  // automatic differentiation does.
834  const unsigned int n_supp = lagrange_support_points.size();
835  const number weight = lagrange_weight;
836  switch (n_derivatives)
837  {
838  default:
839  values[0] = 1.;
840  for (unsigned int d = 1; d <= n_derivatives; ++d)
841  values[d] = 0.;
842  for (unsigned int i = 0; i < n_supp; ++i)
843  {
844  const Number2 v = x - lagrange_support_points[i];
845 
846  // multiply by (x-x_i) and compute action on all derivatives,
847  // too (inspired from automatic differentiation: implement the
848  // product rule for the old value and the new variable 'v',
849  // i.e., expand value v and derivative one). since we reuse a
850  // value from the next lower derivative from the steps before,
851  // need to start from the highest derivative
852  for (unsigned int k = n_derivatives; k > 0; --k)
853  values[k] = (values[k] * v + values[k - 1]);
854  values[0] *= v;
855  }
856  // finally, multiply by the weight in the Lagrange
857  // denominator. Could be done instead of setting values[0] = 1
858  // above, but that gives different accumulation of round-off
859  // errors (multiplication is not associative) compared to when we
860  // computed the weight, and hence a basis function might not be
861  // exactly one at the center point, which is nice to have. We also
862  // multiply derivatives by k! to transform the product p_n =
863  // p^(n)(x)/k! into the actual form of the derivative
864  {
865  number k_factorial = 1;
866  for (unsigned int k = 0; k <= n_derivatives; ++k)
867  {
868  values[k] *= k_factorial * weight;
869  k_factorial *= static_cast<number>(k + 1);
870  }
871  }
872  break;
873 
874  // manually implement case 0 (values only), case 1 (value + first
875  // derivative), and case 2 (up to second derivative) since they
876  // might be called often. then, we can unroll the inner loop and
877  // keep the temporary results as local variables to help the
878  // compiler with the pointer aliasing analysis.
879  case 0:
880  {
881  Number2 value = 1.;
882  for (unsigned int i = 0; i < n_supp; ++i)
883  {
884  const Number2 v = x - lagrange_support_points[i];
885  value *= v;
886  }
887  values[0] = weight * value;
888  break;
889  }
890 
891  case 1:
892  {
893  Number2 value = 1.;
894  Number2 derivative = 0.;
895  for (unsigned int i = 0; i < n_supp; ++i)
896  {
897  const Number2 v = x - lagrange_support_points[i];
898  derivative = derivative * v + value;
899  value *= v;
900  }
901  values[0] = weight * value;
902  values[1] = weight * derivative;
903  break;
904  }
905 
906  case 2:
907  {
908  Number2 value = 1.;
909  Number2 derivative = 0.;
910  Number2 second = 0.;
911  for (unsigned int i = 0; i < n_supp; ++i)
912  {
913  const Number2 v = x - lagrange_support_points[i];
914  second = second * v + derivative;
915  derivative = derivative * v + value;
916  value *= v;
917  }
918  values[0] = weight * value;
919  values[1] = weight * derivative;
920  values[2] = static_cast<number>(2) * weight * second;
921  break;
922  }
923  }
924  return;
925  }
926 
927  Assert(coefficients.size() > 0, ExcEmptyObject());
928 
929  // if derivatives are needed, then do it properly by the full
930  // Horner scheme
931  const unsigned int m = coefficients.size();
932  std::vector<Number2> a(coefficients.size());
933  std::copy(coefficients.begin(), coefficients.end(), a.begin());
934  unsigned int j_factorial = 1;
935 
936  // loop over all requested derivatives. note that derivatives @p{j>m} are
937  // necessarily zero, as they differentiate the polynomial more often than
938  // the highest power is
939  const unsigned int min_valuessize_m = std::min(n_derivatives + 1, m);
940  for (unsigned int j = 0; j < min_valuessize_m; ++j)
941  {
942  for (int k = m - 2; k >= static_cast<int>(j); --k)
943  a[k] += x * a[k + 1];
944  values[j] = static_cast<number>(j_factorial) * a[j];
945 
946  j_factorial *= j + 1;
947  }
948 
949  // fill higher derivatives by zero
950  for (unsigned int j = min_valuessize_m; j <= n_derivatives; ++j)
951  values[j] = 0.;
952  }
953 
954 
955 
956  template <typename number>
957  template <class Archive>
958  inline void
959  Polynomial<number>::serialize(Archive &ar, const unsigned int)
960  {
961  // forward to serialization function in the base class.
962  ar &static_cast<Subscriptor &>(*this);
963  ar &coefficients;
964  ar &in_lagrange_product_form;
965  ar &lagrange_support_points;
966  ar &lagrange_weight;
967  }
968 
969 
970 
971  template <typename Number>
972  Number
973  jacobi_polynomial_value(const unsigned int degree,
974  const int alpha,
975  const int beta,
976  const Number x)
977  {
978  Assert(alpha >= 0 && beta >= 0,
979  ExcNotImplemented("Negative alpha/beta coefficients not supported"));
980  // the Jacobi polynomial is evaluated using a recursion formula.
981  Number p0, p1;
982 
983  // The recursion formula is defined for the interval [-1, 1], so rescale
984  // to that interval here
985  const Number xeval = Number(-1) + 2. * x;
986 
987  // initial values P_0(x), P_1(x):
988  p0 = 1.0;
989  if (degree == 0)
990  return p0;
991  p1 = ((alpha + beta + 2) * xeval + (alpha - beta)) / 2;
992  if (degree == 1)
993  return p1;
994 
995  for (unsigned int i = 1; i < degree; ++i)
996  {
997  const Number v = 2 * i + (alpha + beta);
998  const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
999  const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
1000  const Number a3 = v * (v + 1) * (v + 2);
1001  const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
1002 
1003  const Number pn = ((a2 + a3 * xeval) * p1 - a4 * p0) / a1;
1004  p0 = p1;
1005  p1 = pn;
1006  }
1007  return p1;
1008  }
1009 
1010 
1011 
1012  template <typename Number>
1013  std::vector<Number>
1014  jacobi_polynomial_roots(const unsigned int degree,
1015  const int alpha,
1016  const int beta)
1017  {
1018  std::vector<Number> x(degree, 0.5);
1019 
1020  // compute zeros with a Newton algorithm.
1021 
1022  // Set tolerance. For long double we might not always get the additional
1023  // precision in a run time environment (e.g. with valgrind), so we must
1024  // limit the tolerance to double. Since we do a Newton iteration, doing
1025  // one more iteration after the residual has indicated convergence will be
1026  // enough for all number types due to the quadratic convergence of
1027  // Newton's method
1028 
1029  const Number tolerance =
1030  4 * std::max(static_cast<Number>(std::numeric_limits<double>::epsilon()),
1032 
1033  // The following implementation follows closely the one given in the
1034  // appendix of the book by Karniadakis and Sherwin: Spectral/hp element
1035  // methods for computational fluid dynamics (Oxford University Press,
1036  // 2005)
1037 
1038  // If symmetric, we only need to compute the half of points
1039  const unsigned int n_points = (alpha == beta ? degree / 2 : degree);
1040  for (unsigned int k = 0; k < n_points; ++k)
1041  {
1042  // we take the zeros of the Chebyshev polynomial (alpha=beta=-0.5) as
1043  // initial values, corrected by the initial value
1044  Number r = 0.5 - 0.5 * std::cos(static_cast<Number>(2 * k + 1) /
1045  (2 * degree) * numbers::PI);
1046  if (k > 0)
1047  r = (r + x[k - 1]) / 2;
1048 
1049  unsigned int converged = numbers::invalid_unsigned_int;
1050  for (unsigned int it = 1; it < 1000; ++it)
1051  {
1052  Number s = 0.;
1053  for (unsigned int i = 0; i < k; ++i)
1054  s += 1. / (r - x[i]);
1055 
1056  // derivative of P_n^{alpha,beta}, rescaled to [0, 1]
1057  const Number J_x =
1058  (alpha + beta + degree + 1) *
1059  jacobi_polynomial_value(degree - 1, alpha + 1, beta + 1, r);
1060 
1061  // value of P_n^{alpha,beta}
1062  const Number f = jacobi_polynomial_value(degree, alpha, beta, r);
1063  const Number delta = f / (f * s - J_x);
1064  r += delta;
1065  if (converged == numbers::invalid_unsigned_int &&
1066  std::abs(delta) < tolerance)
1067  converged = it;
1068 
1069  // do one more iteration to ensure accuracy also for tighter
1070  // types than double (e.g. long double)
1071  if (it == converged + 1)
1072  break;
1073  }
1074 
1076  ExcMessage("Newton iteration for zero of Jacobi polynomial "
1077  "did not converge."));
1078 
1079  x[k] = r;
1080  }
1081 
1082  // in case we assumed symmetry, fill up the missing values
1083  for (unsigned int k = n_points; k < degree; ++k)
1084  x[k] = 1.0 - x[degree - k - 1];
1085 
1086  return x;
1087  }
1088 
1089 } // namespace Polynomials
1091 
1092 #endif
static const unsigned int invalid_unsigned_int
Definition: types.h:196
void serialize(Archive &ar, const unsigned int version)
Definition: polynomial.h:959
Polynomial< number > & operator*=(const double s)
Definition: polynomial.cc:203
void scale(const number factor)
Definition: polynomial.cc:165
void transform_into_standard_form()
Definition: polynomial.cc:111
bool operator==(const Polynomial< number > &p) const
Definition: polynomial.cc:346
Polynomial< number > derivative() const
Definition: polynomial.cc:458
void shift(const number2 offset)
Definition: polynomial.cc:439
unsigned int degree() const
Definition: polynomial.h:777
Polynomial< number > & operator-=(const Polynomial< number > &p)
Definition: polynomial.cc:310
static std::vector< std::unique_ptr< const std::vector< double > > > recursive_coefficients
Definition: polynomial.h:551
Point< 2 > second
Definition: grid_out.cc:4356
std::vector< Number > jacobi_polynomial_roots(const unsigned int degree, const int alpha, const int beta)
Definition: polynomial.h:1014
number value(const number x) const
Definition: polynomial.h:794
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
void print(std::ostream &out) const
Definition: polynomial.cc:514
Number jacobi_polynomial_value(const unsigned int degree, const int alpha, const int beta, const Number x)
Definition: polynomial.h:973
static ::ExceptionBase & ExcMessage(std::string arg1)
static void multiply(std::vector< number > &coefficients, const number factor)
Definition: polynomial.cc:190
#define Assert(cond, exc)
Definition: exceptions.h:1466
Polynomial< number > & operator+=(const Polynomial< number > &p)
Definition: polynomial.cc:268
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:380
virtual std::size_t memory_consumption() const
Definition: polynomial.cc:533
std::vector< number > coefficients
Definition: polynomial.h:279
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Polynomial< number > primitive() const
Definition: polynomial.cc:487
std::vector< number > lagrange_support_points
Definition: polynomial.h:291
static constexpr double PI
Definition: numbers.h:231
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:379
T min(const T &t, const MPI_Comm &mpi_communicator)
static ::ExceptionBase & ExcEmptyObject()
static ::ExceptionBase & ExcNotImplemented()
void copy(const T *begin, const T *end, U *dest)
T max(const T &t, const MPI_Comm &mpi_communicator)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
Definition: polynomial.cc:701