Reference documentation for deal.II version Git 1ecc23629d 2021-05-18 09:57:04 +0200
\(\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 <memory>
28 #include <vector>
29 
31 
41 namespace Polynomials
42 {
62  template <typename number>
63  class Polynomial : public Subscriptor
64  {
65  public:
74  Polynomial(const std::vector<number> &coefficients);
75 
79  Polynomial(const unsigned int n);
80 
88  Polynomial(const std::vector<Point<1>> &lagrange_support_points,
89  const unsigned int evaluation_point);
90 
94  Polynomial();
95 
105  number
106  value(const number x) const;
107 
118  void
119  value(const number x, std::vector<number> &values) const;
120 
139  template <typename Number2>
140  void
141  value(const Number2 x,
142  const unsigned int n_derivatives,
143  Number2 * values) const;
144 
150  unsigned int
151  degree() const;
152 
160  void
161  scale(const number factor);
162 
178  template <typename number2>
179  void
180  shift(const number2 offset);
181 
186  derivative() const;
187 
193  primitive() const;
194 
199  operator*=(const double s);
200 
205  operator*=(const Polynomial<number> &p);
206 
211  operator+=(const Polynomial<number> &p);
212 
217  operator-=(const Polynomial<number> &p);
218 
222  bool
223  operator==(const Polynomial<number> &p) const;
224 
228  void
229  print(std::ostream &out) const;
230 
236  template <class Archive>
237  void
238  serialize(Archive &ar, const unsigned int version);
239 
243  virtual std::size_t
244  memory_consumption() const;
245 
246  protected:
250  static void
251  scale(std::vector<number> &coefficients, const number factor);
252 
256  template <typename number2>
257  static void
258  shift(std::vector<number> &coefficients, const number2 shift);
259 
263  static void
264  multiply(std::vector<number> &coefficients, const number factor);
265 
271  void
273 
282  std::vector<number> coefficients;
283 
289 
294  std::vector<number> lagrange_support_points;
295 
301  };
302 
303 
308  template <typename number>
309  class Monomial : public Polynomial<number>
310  {
311  public:
316  Monomial(const unsigned int n, const double coefficient = 1.);
317 
324  static std::vector<Polynomial<number>>
325  generate_complete_basis(const unsigned int degree);
326 
327  private:
331  static std::vector<number>
332  make_vector(unsigned int n, const double coefficient);
333  };
334 
335 
351  class LagrangeEquidistant : public Polynomial<double>
352  {
353  public:
359  LagrangeEquidistant(const unsigned int n, const unsigned int support_point);
360 
369  static std::vector<Polynomial<double>>
370  generate_complete_basis(const unsigned int degree);
371 
372  private:
377  static void
378  compute_coefficients(const unsigned int n,
379  const unsigned int support_point,
380  std::vector<double> &a);
381  };
382 
383 
384 
391  std::vector<Polynomial<double>>
392  generate_complete_Lagrange_basis(const std::vector<Point<1>> &points);
393 
394 
395 
409  class Legendre : public Polynomial<double>
410  {
411  public:
415  Legendre(const unsigned int p);
416 
423  static std::vector<Polynomial<double>>
424  generate_complete_basis(const unsigned int degree);
425  };
426 
446  class Lobatto : public Polynomial<double>
447  {
448  public:
453  Lobatto(const unsigned int p = 0);
454 
459  static std::vector<Polynomial<double>>
460  generate_complete_basis(const unsigned int p);
461 
462  private:
466  std::vector<double>
467  compute_coefficients(const unsigned int p);
468  };
469 
470 
471 
509  class Hierarchical : public Polynomial<double>
510  {
511  public:
516  Hierarchical(const unsigned int p);
517 
528  static std::vector<Polynomial<double>>
529  generate_complete_basis(const unsigned int degree);
530 
531  private:
535  static void
536  compute_coefficients(const unsigned int p);
537 
542  static const std::vector<double> &
543  get_coefficients(const unsigned int p);
544 
553  static std::vector<std::unique_ptr<const std::vector<double>>>
555  };
556 
557 
558 
586  class HermiteInterpolation : public Polynomial<double>
587  {
588  public:
593  HermiteInterpolation(const unsigned int p);
594 
600  static std::vector<Polynomial<double>>
601  generate_complete_basis(const unsigned int p);
602  };
603 
604 
605 
707  class HermiteLikeInterpolation : public Polynomial<double>
708  {
709  public:
714  HermiteLikeInterpolation(const unsigned int degree,
715  const unsigned int index);
716 
721  static std::vector<Polynomial<double>>
722  generate_complete_basis(const unsigned int degree);
723  };
724 
725 
726 
727  /*
728  * Evaluate a Jacobi polynomial @f$ P_n^{\alpha, \beta}(x) @f$ specified by the
729  * parameters @p alpha, @p beta, @p n, where @p n is the degree of the
730  * Jacobi polynomial.
731  *
732  * @note The Jacobi polynomials are not orthonormal and are defined on the
733  * unit interval @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often
734  * used in literature. @p x is the point of evaluation.
735  */
736  template <typename Number>
737  Number
738  jacobi_polynomial_value(const unsigned int degree,
739  const int alpha,
740  const int beta,
741  const Number x);
742 
743 
756  template <typename Number>
757  std::vector<Number>
758  jacobi_polynomial_roots(const unsigned int degree,
759  const int alpha,
760  const int beta);
761 } // namespace Polynomials
762 
763 
766 /* -------------------------- inline functions --------------------- */
767 
768 namespace Polynomials
769 {
770  template <typename number>
772  : in_lagrange_product_form(false)
773  , lagrange_weight(1.)
774  {}
775 
776 
777 
778  template <typename number>
779  inline unsigned int
781  {
782  if (in_lagrange_product_form == true)
783  {
784  return lagrange_support_points.size();
785  }
786  else
787  {
788  Assert(coefficients.size() > 0, ExcEmptyObject());
789  return coefficients.size() - 1;
790  }
791  }
792 
793 
794 
795  template <typename number>
796  inline number
797  Polynomial<number>::value(const number x) const
798  {
799  if (in_lagrange_product_form == false)
800  {
801  Assert(coefficients.size() > 0, ExcEmptyObject());
802 
803  // Horner scheme
804  const unsigned int m = coefficients.size();
805  number value = coefficients.back();
806  for (int k = m - 2; k >= 0; --k)
807  value = value * x + coefficients[k];
808  return value;
809  }
810  else
811  {
812  // direct evaluation of Lagrange polynomial
813  const unsigned int m = lagrange_support_points.size();
814  number value = 1.;
815  for (unsigned int j = 0; j < m; ++j)
816  value *= x - lagrange_support_points[j];
817  value *= lagrange_weight;
818  return value;
819  }
820  }
821 
822 
823 
824  template <typename number>
825  template <typename Number2>
826  inline void
827  Polynomial<number>::value(const Number2 x,
828  const unsigned int n_derivatives,
829  Number2 * values) const
830  {
831  // evaluate Lagrange polynomial and derivatives
832  if (in_lagrange_product_form == true)
833  {
834  // to compute the value and all derivatives of a polynomial of the
835  // form (x-x_1)*(x-x_2)*...*(x-x_n), expand the derivatives like
836  // automatic differentiation does.
837  const unsigned int n_supp = lagrange_support_points.size();
838  const number weight = lagrange_weight;
839  switch (n_derivatives)
840  {
841  default:
842  values[0] = 1.;
843  for (unsigned int d = 1; d <= n_derivatives; ++d)
844  values[d] = 0.;
845  for (unsigned int i = 0; i < n_supp; ++i)
846  {
847  const Number2 v = x - lagrange_support_points[i];
848 
849  // multiply by (x-x_i) and compute action on all derivatives,
850  // too (inspired from automatic differentiation: implement the
851  // product rule for the old value and the new variable 'v',
852  // i.e., expand value v and derivative one). since we reuse a
853  // value from the next lower derivative from the steps before,
854  // need to start from the highest derivative
855  for (unsigned int k = n_derivatives; k > 0; --k)
856  values[k] = (values[k] * v + values[k - 1]);
857  values[0] *= v;
858  }
859  // finally, multiply by the weight in the Lagrange
860  // denominator. Could be done instead of setting values[0] = 1
861  // above, but that gives different accumulation of round-off
862  // errors (multiplication is not associative) compared to when we
863  // computed the weight, and hence a basis function might not be
864  // exactly one at the center point, which is nice to have. We also
865  // multiply derivatives by k! to transform the product p_n =
866  // p^(n)(x)/k! into the actual form of the derivative
867  {
868  number k_factorial = 1;
869  for (unsigned int k = 0; k <= n_derivatives; ++k)
870  {
871  values[k] *= k_factorial * weight;
872  k_factorial *= static_cast<number>(k + 1);
873  }
874  }
875  break;
876 
877  // manually implement case 0 (values only), case 1 (value + first
878  // derivative), and case 2 (up to second derivative) since they
879  // might be called often. then, we can unroll the inner loop and
880  // keep the temporary results as local variables to help the
881  // compiler with the pointer aliasing analysis.
882  case 0:
883  {
884  Number2 value = 1.;
885  for (unsigned int i = 0; i < n_supp; ++i)
886  {
887  const Number2 v = x - lagrange_support_points[i];
888  value *= v;
889  }
890  values[0] = weight * value;
891  break;
892  }
893 
894  case 1:
895  {
896  Number2 value = 1.;
897  Number2 derivative = 0.;
898  for (unsigned int i = 0; i < n_supp; ++i)
899  {
900  const Number2 v = x - lagrange_support_points[i];
901  derivative = derivative * v + value;
902  value *= v;
903  }
904  values[0] = weight * value;
905  values[1] = weight * derivative;
906  break;
907  }
908 
909  case 2:
910  {
911  Number2 value = 1.;
912  Number2 derivative = 0.;
913  Number2 second = 0.;
914  for (unsigned int i = 0; i < n_supp; ++i)
915  {
916  const Number2 v = x - lagrange_support_points[i];
917  second = second * v + derivative;
918  derivative = derivative * v + value;
919  value *= v;
920  }
921  values[0] = weight * value;
922  values[1] = weight * derivative;
923  values[2] = static_cast<number>(2) * weight * second;
924  break;
925  }
926  }
927  return;
928  }
929 
930  Assert(coefficients.size() > 0, ExcEmptyObject());
931 
932  // if derivatives are needed, then do it properly by the full
933  // Horner scheme
934  const unsigned int m = coefficients.size();
935  std::vector<Number2> a(coefficients.size());
936  std::copy(coefficients.begin(), coefficients.end(), a.begin());
937  unsigned int j_factorial = 1;
938 
939  // loop over all requested derivatives. note that derivatives @p{j>m} are
940  // necessarily zero, as they differentiate the polynomial more often than
941  // the highest power is
942  const unsigned int min_valuessize_m = std::min(n_derivatives + 1, m);
943  for (unsigned int j = 0; j < min_valuessize_m; ++j)
944  {
945  for (int k = m - 2; k >= static_cast<int>(j); --k)
946  a[k] += x * a[k + 1];
947  values[j] = static_cast<number>(j_factorial) * a[j];
948 
949  j_factorial *= j + 1;
950  }
951 
952  // fill higher derivatives by zero
953  for (unsigned int j = min_valuessize_m; j <= n_derivatives; ++j)
954  values[j] = 0.;
955  }
956 
957 
958 
959  template <typename number>
960  template <class Archive>
961  inline void
962  Polynomial<number>::serialize(Archive &ar, const unsigned int)
963  {
964  // forward to serialization function in the base class.
965  ar &static_cast<Subscriptor &>(*this);
966  ar &coefficients;
967  ar &in_lagrange_product_form;
968  ar &lagrange_support_points;
969  ar &lagrange_weight;
970  }
971 
972 
973 
974  template <typename Number>
975  Number
976  jacobi_polynomial_value(const unsigned int degree,
977  const int alpha,
978  const int beta,
979  const Number x)
980  {
981  Assert(alpha >= 0 && beta >= 0,
982  ExcNotImplemented("Negative alpha/beta coefficients not supported"));
983  // the Jacobi polynomial is evaluated using a recursion formula.
984  Number p0, p1;
985 
986  // The recursion formula is defined for the interval [-1, 1], so rescale
987  // to that interval here
988  const Number xeval = Number(-1) + 2. * x;
989 
990  // initial values P_0(x), P_1(x):
991  p0 = 1.0;
992  if (degree == 0)
993  return p0;
994  p1 = ((alpha + beta + 2) * xeval + (alpha - beta)) / 2;
995  if (degree == 1)
996  return p1;
997 
998  for (unsigned int i = 1; i < degree; ++i)
999  {
1000  const Number v = 2 * i + (alpha + beta);
1001  const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
1002  const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
1003  const Number a3 = v * (v + 1) * (v + 2);
1004  const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
1005 
1006  const Number pn = ((a2 + a3 * xeval) * p1 - a4 * p0) / a1;
1007  p0 = p1;
1008  p1 = pn;
1009  }
1010  return p1;
1011  }
1012 
1013 
1014 
1015  template <typename Number>
1016  std::vector<Number>
1017  jacobi_polynomial_roots(const unsigned int degree,
1018  const int alpha,
1019  const int beta)
1020  {
1021  std::vector<Number> x(degree, 0.5);
1022 
1023  // compute zeros with a Newton algorithm.
1024 
1025  // Set tolerance. For long double we might not always get the additional
1026  // precision in a run time environment (e.g. with valgrind), so we must
1027  // limit the tolerance to double. Since we do a Newton iteration, doing
1028  // one more iteration after the residual has indicated convergence will be
1029  // enough for all number types due to the quadratic convergence of
1030  // Newton's method
1031 
1032  const Number tolerance =
1033  4 * std::max(static_cast<Number>(std::numeric_limits<double>::epsilon()),
1035 
1036  // The following implementation follows closely the one given in the
1037  // appendix of the book by Karniadakis and Sherwin: Spectral/hp element
1038  // methods for computational fluid dynamics (Oxford University Press,
1039  // 2005)
1040 
1041  // If symmetric, we only need to compute the half of points
1042  const unsigned int n_points = (alpha == beta ? degree / 2 : degree);
1043  for (unsigned int k = 0; k < n_points; ++k)
1044  {
1045  // we take the zeros of the Chebyshev polynomial (alpha=beta=-0.5) as
1046  // initial values, corrected by the initial value
1047  Number r = 0.5 - 0.5 * std::cos(static_cast<Number>(2 * k + 1) /
1048  (2 * degree) * numbers::PI);
1049  if (k > 0)
1050  r = (r + x[k - 1]) / 2;
1051 
1052  unsigned int converged = numbers::invalid_unsigned_int;
1053  for (unsigned int it = 1; it < 1000; ++it)
1054  {
1055  Number s = 0.;
1056  for (unsigned int i = 0; i < k; ++i)
1057  s += 1. / (r - x[i]);
1058 
1059  // derivative of P_n^{alpha,beta}, rescaled to [0, 1]
1060  const Number J_x =
1061  (alpha + beta + degree + 1) *
1062  jacobi_polynomial_value(degree - 1, alpha + 1, beta + 1, r);
1063 
1064  // value of P_n^{alpha,beta}
1065  const Number f = jacobi_polynomial_value(degree, alpha, beta, r);
1066  const Number delta = f / (f * s - J_x);
1067  r += delta;
1068  if (converged == numbers::invalid_unsigned_int &&
1069  std::abs(delta) < tolerance)
1070  converged = it;
1071 
1072  // do one more iteration to ensure accuracy also for tighter
1073  // types than double (e.g. long double)
1074  if (it == converged + 1)
1075  break;
1076  }
1077 
1079  ExcMessage("Newton iteration for zero of Jacobi polynomial "
1080  "did not converge."));
1081 
1082  x[k] = r;
1083  }
1084 
1085  // in case we assumed symmetry, fill up the missing values
1086  for (unsigned int k = n_points; k < degree; ++k)
1087  x[k] = 1.0 - x[degree - k - 1];
1088 
1089  return x;
1090  }
1091 
1092 } // namespace Polynomials
1094 
1095 #endif
static const unsigned int invalid_unsigned_int
Definition: types.h:196
void serialize(Archive &ar, const unsigned int version)
Definition: polynomial.h:962
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:780
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:554
Point< 2 > second
Definition: grid_out.cc:4588
std::vector< Number > jacobi_polynomial_roots(const unsigned int degree, const int alpha, const int beta)
Definition: polynomial.h:1017
number value(const number x) const
Definition: polynomial.h:797
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:976
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:1465
Polynomial< number > & operator+=(const Polynomial< number > &p)
Definition: polynomial.cc:268
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
virtual std::size_t memory_consumption() const
Definition: polynomial.cc:533
std::vector< number > coefficients
Definition: polynomial.h:282
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:294
static constexpr double PI
Definition: numbers.h:231
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
static ::ExceptionBase & ExcEmptyObject()
static ::ExceptionBase & ExcNotImplemented()
void copy(const T *begin, const T *end, U *dest)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
Definition: polynomial.cc:701