Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
polynomial.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2018 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 
23 #include <deal.II/base/exceptions.h>
24 #include <deal.II/base/point.h>
25 #include <deal.II/base/subscriptor.h>
26 
27 #include <memory>
28 #include <vector>
29 
30 DEAL_II_NAMESPACE_OPEN
31 
41 namespace Polynomials
42 {
61  template <typename number>
62  class Polynomial : public Subscriptor
63  {
64  public:
73  Polynomial(const std::vector<number> &coefficients);
74 
78  Polynomial(const unsigned int n);
79 
87  Polynomial(const std::vector<Point<1>> &lagrange_support_points,
88  const unsigned int evaluation_point);
89 
93  Polynomial();
94 
102  number
103  value(const number x) const;
104 
115  void
116  value(const number x, std::vector<number> &values) const;
117 
129  void
130  value(const number x,
131  const unsigned int n_derivatives,
132  number * values) const;
133 
139  unsigned int
140  degree() const;
141 
149  void
150  scale(const number factor);
151 
167  template <typename number2>
168  void
169  shift(const number2 offset);
170 
175  derivative() const;
176 
182  primitive() const;
183 
188  operator*=(const double s);
189 
194  operator*=(const Polynomial<number> &p);
195 
200  operator+=(const Polynomial<number> &p);
201 
206  operator-=(const Polynomial<number> &p);
207 
211  bool
212  operator==(const Polynomial<number> &p) const;
213 
217  void
218  print(std::ostream &out) const;
219 
224  template <class Archive>
225  void
226  serialize(Archive &ar, const unsigned int version);
227 
228  protected:
232  static void
233  scale(std::vector<number> &coefficients, const number factor);
234 
238  template <typename number2>
239  static void
240  shift(std::vector<number> &coefficients, const number2 shift);
241 
245  static void
246  multiply(std::vector<number> &coefficients, const number factor);
247 
253  void
255 
264  std::vector<number> coefficients;
265 
271 
276  std::vector<number> lagrange_support_points;
277 
283  };
284 
285 
292  template <typename number>
293  class Monomial : public Polynomial<number>
294  {
295  public:
300  Monomial(const unsigned int n, const double coefficient = 1.);
301 
308  static std::vector<Polynomial<number>>
309  generate_complete_basis(const unsigned int degree);
310 
311  private:
315  static std::vector<number>
316  make_vector(unsigned int n, const double coefficient);
317  };
318 
319 
337  class LagrangeEquidistant : public Polynomial<double>
338  {
339  public:
345  LagrangeEquidistant(const unsigned int n, const unsigned int support_point);
346 
355  static std::vector<Polynomial<double>>
356  generate_complete_basis(const unsigned int degree);
357 
358  private:
363  static void
364  compute_coefficients(const unsigned int n,
365  const unsigned int support_point,
366  std::vector<double> &a);
367  };
368 
369 
370 
377  std::vector<Polynomial<double>>
378  generate_complete_Lagrange_basis(const std::vector<Point<1>> &points);
379 
380 
381 
397  class Legendre : public Polynomial<double>
398  {
399  public:
403  Legendre(const unsigned int p);
404 
411  static std::vector<Polynomial<double>>
412  generate_complete_basis(const unsigned int degree);
413  };
414 
436  class Lobatto : public Polynomial<double>
437  {
438  public:
443  Lobatto(const unsigned int p = 0);
444 
449  static std::vector<Polynomial<double>>
450  generate_complete_basis(const unsigned int p);
451 
452  private:
456  std::vector<double>
457  compute_coefficients(const unsigned int p);
458  };
459 
460 
461 
501  class Hierarchical : public Polynomial<double>
502  {
503  public:
508  Hierarchical(const unsigned int p);
509 
520  static std::vector<Polynomial<double>>
521  generate_complete_basis(const unsigned int degree);
522 
523  private:
527  static void
528  compute_coefficients(const unsigned int p);
529 
534  static const std::vector<double> &
535  get_coefficients(const unsigned int p);
536 
545  static std::vector<std::unique_ptr<const std::vector<double>>>
547  };
548 
549 
550 
581  class HermiteInterpolation : public Polynomial<double>
582  {
583  public:
588  HermiteInterpolation(const unsigned int p);
589 
595  static std::vector<Polynomial<double>>
596  generate_complete_basis(const unsigned int p);
597  };
598 
599 
600 
702  class HermiteLikeInterpolation : public Polynomial<double>
703  {
704  public:
709  HermiteLikeInterpolation(const unsigned int degree,
710  const unsigned int index);
711 
716  static std::vector<Polynomial<double>>
717  generate_complete_basis(const unsigned int degree);
718  };
719 
720 
721 
722  /*
723  * Evaluate a Jacobi polynomial @f$ P_n^{\alpha, \beta}(x) @f$ specified by the
724  * parameters @p alpha, @p beta, @p n, where @p n is the degree of the
725  * Jacobi polynomial.
726  *
727  * @note The Jacobi polynomials are not orthonormal and are defined on the
728  * unit interval @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often
729  * used in literature. @p x is the point of evaluation.
730  */
731  template <typename Number>
732  Number
733  jacobi_polynomial_value(const unsigned int degree,
734  const int alpha,
735  const int beta,
736  const Number x);
737 
738 
751  template <typename Number>
752  std::vector<Number>
753  jacobi_polynomial_roots(const unsigned int degree,
754  const int alpha,
755  const int beta);
756 } // namespace Polynomials
757 
758 
761 /* -------------------------- inline functions --------------------- */
762 
763 namespace Polynomials
764 {
765  template <typename number>
767  : in_lagrange_product_form(false)
768  , lagrange_weight(1.)
769  {}
770 
771 
772 
773  template <typename number>
774  inline unsigned int
776  {
777  if (in_lagrange_product_form == true)
778  {
779  return lagrange_support_points.size();
780  }
781  else
782  {
783  Assert(coefficients.size() > 0, ExcEmptyObject());
784  return coefficients.size() - 1;
785  }
786  }
787 
788 
789 
790  template <typename number>
791  inline number
792  Polynomial<number>::value(const number x) const
793  {
794  if (in_lagrange_product_form == false)
795  {
796  Assert(coefficients.size() > 0, ExcEmptyObject());
797 
798  // Horner scheme
799  const unsigned int m = coefficients.size();
800  number value = coefficients.back();
801  for (int k = m - 2; k >= 0; --k)
802  value = value * x + coefficients[k];
803  return value;
804  }
805  else
806  {
807  // direct evaluation of Lagrange polynomial
808  const unsigned int m = lagrange_support_points.size();
809  number value = 1.;
810  for (unsigned int j = 0; j < m; ++j)
811  value *= x - lagrange_support_points[j];
812  value *= lagrange_weight;
813  return value;
814  }
815  }
816 
817 
818 
819  template <typename number>
820  template <class Archive>
821  inline void
822  Polynomial<number>::serialize(Archive &ar, const unsigned int)
823  {
824  // forward to serialization function in the base class.
825  ar &static_cast<Subscriptor &>(*this);
826  ar &coefficients;
827  ar &in_lagrange_product_form;
828  ar &lagrange_support_points;
829  ar &lagrange_weight;
830  }
831 
832 
833 
834  template <typename Number>
835  Number
836  jacobi_polynomial_value(const unsigned int degree,
837  const int alpha,
838  const int beta,
839  const Number x)
840  {
841  Assert(alpha >= 0 && beta >= 0,
842  ExcNotImplemented("Negative alpha/beta coefficients not supported"));
843  // the Jacobi polynomial is evaluated using a recursion formula.
844  Number p0, p1;
845 
846  // The recursion formula is defined for the interval [-1, 1], so rescale
847  // to that interval here
848  const Number xeval = Number(-1) + 2. * x;
849 
850  // initial values P_0(x), P_1(x):
851  p0 = 1.0;
852  if (degree == 0)
853  return p0;
854  p1 = ((alpha + beta + 2) * xeval + (alpha - beta)) / 2;
855  if (degree == 1)
856  return p1;
857 
858  for (unsigned int i = 1; i < degree; ++i)
859  {
860  const Number v = 2 * i + (alpha + beta);
861  const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
862  const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
863  const Number a3 = v * (v + 1) * (v + 2);
864  const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
865 
866  const Number pn = ((a2 + a3 * xeval) * p1 - a4 * p0) / a1;
867  p0 = p1;
868  p1 = pn;
869  }
870  return p1;
871  }
872 
873 
874 
875  template <typename Number>
876  std::vector<Number>
877  jacobi_polynomial_roots(const unsigned int degree,
878  const int alpha,
879  const int beta)
880  {
881  std::vector<Number> x(degree, 0.5);
882 
883  // compute zeros with a Newton algorithm.
884 
885  // Set tolerance. For long double we might not always get the additional
886  // precision in a run time environment (e.g. with valgrind), so we must
887  // limit the tolerance to double. Since we do a Newton iteration, doing
888  // one more iteration after the residual has indicated convergence will be
889  // enough for all number types due to the quadratic convergence of
890  // Newton's method
891 
892  const Number tolerance =
893  4 * std::max(static_cast<Number>(std::numeric_limits<double>::epsilon()),
894  std::numeric_limits<Number>::epsilon());
895 
896  // The following implementation follows closely the one given in the
897  // appendix of the book by Karniadakis and Sherwin: Spectral/hp element
898  // methods for computational fluid dynamics (Oxford University Press,
899  // 2005)
900 
901  // If symmetric, we only need to compute the half of points
902  const unsigned int n_points = (alpha == beta ? degree / 2 : degree);
903  for (unsigned int k = 0; k < n_points; ++k)
904  {
905  // we take the zeros of the Chebyshev polynomial (alpha=beta=-0.5) as
906  // initial values, corrected by the initial value
907  Number r = 0.5 - 0.5 * std::cos(static_cast<Number>(2 * k + 1) /
908  (2 * degree) * numbers::PI);
909  if (k > 0)
910  r = (r + x[k - 1]) / 2;
911 
912  unsigned int converged = numbers::invalid_unsigned_int;
913  for (unsigned int it = 1; it < 1000; ++it)
914  {
915  Number s = 0.;
916  for (unsigned int i = 0; i < k; ++i)
917  s += 1. / (r - x[i]);
918 
919  // derivative of P_n^{alpha,beta}, rescaled to [0, 1]
920  const Number J_x =
921  (alpha + beta + degree + 1) *
922  jacobi_polynomial_value(degree - 1, alpha + 1, beta + 1, r);
923 
924  // value of P_n^{alpha,beta}
925  const Number f = jacobi_polynomial_value(degree, alpha, beta, r);
926  const Number delta = f / (f * s - J_x);
927  r += delta;
928  if (converged == numbers::invalid_unsigned_int &&
929  std::abs(delta) < tolerance)
930  converged = it;
931 
932  // do one more iteration to ensure accuracy also for tighter
933  // types than double (e.g. long double)
934  if (it == converged + 1)
935  break;
936  }
937 
939  ExcMessage("Newton iteration for zero of Jacobi polynomial "
940  "did not converge."));
941 
942  x[k] = r;
943  }
944 
945  // in case we assumed symmetry, fill up the missing values
946  for (unsigned int k = n_points; k < degree; ++k)
947  x[k] = 1.0 - x[degree - k - 1];
948 
949  return x;
950  }
951 
952 } // namespace Polynomials
953 DEAL_II_NAMESPACE_CLOSE
954 
955 #endif
static const unsigned int invalid_unsigned_int
Definition: types.h:173
void serialize(Archive &ar, const unsigned int version)
Definition: polynomial.h:822
Polynomial< number > & operator*=(const double s)
Definition: polynomial.cc:335
void scale(const number factor)
Definition: polynomial.cc:297
void transform_into_standard_form()
Definition: polynomial.cc:243
bool operator==(const Polynomial< number > &p) const
Definition: polynomial.cc:478
LagrangeEquidistant(const unsigned int n, const unsigned int support_point)
Definition: polynomial.cc:717
Polynomial< number > derivative() const
Definition: polynomial.cc:590
void shift(const number2 offset)
Definition: polynomial.cc:571
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
Definition: polynomial.cc:971
unsigned int degree() const
Definition: polynomial.h:775
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:1150
Legendre(const unsigned int p)
Definition: polynomial.cc:839
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:866
Polynomial< number > & operator-=(const Polynomial< number > &p)
Definition: polynomial.cc:442
static std::vector< std::unique_ptr< const std::vector< double > > > recursive_coefficients
Definition: polynomial.h:546
std::vector< Number > jacobi_polynomial_roots(const unsigned int degree, const int alpha, const int beta)
Definition: polynomial.h:877
static void compute_coefficients(const unsigned int n, const unsigned int support_point, std::vector< double > &a)
Definition: polynomial.cc:741
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:800
Monomial(const unsigned int n, const double coefficient=1.)
Definition: polynomial.cc:678
number value(const number x) const
Definition: polynomial.h:792
Definition: point.h:110
void print(std::ostream &out) const
Definition: polynomial.cc:646
std::vector< double > compute_coefficients(const unsigned int p)
Definition: polynomial.cc:885
static ::ExceptionBase & ExcMessage(std::string arg1)
static void multiply(std::vector< number > &coefficients, const number factor)
Definition: polynomial.cc:322
#define Assert(cond, exc)
Definition: exceptions.h:1407
Polynomial< number > & operator+=(const Polynomial< number > &p)
Definition: polynomial.cc:400
HermiteInterpolation(const unsigned int p)
Definition: polynomial.cc:1178
Lobatto(const unsigned int p=0)
Definition: polynomial.cc:880
static std::vector< Polynomial< number > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:686
static std::vector< number > make_vector(unsigned int n, const double coefficient)
Definition: polynomial.cc:668
Hierarchical(const unsigned int p)
Definition: polynomial.cc:993
HermiteLikeInterpolation(const unsigned int degree, const unsigned int index)
Definition: polynomial.cc:1296
std::vector< number > coefficients
Definition: polynomial.h:264
Polynomial< number > primitive() const
Definition: polynomial.cc:619
static void compute_coefficients(const unsigned int p)
Definition: polynomial.cc:1000
static const std::vector< double > & get_coefficients(const unsigned int p)
Definition: polynomial.cc:1134
std::vector< number > lagrange_support_points
Definition: polynomial.h:276
static constexpr double PI
Definition: numbers.h:146
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:1524
static ::ExceptionBase & ExcEmptyObject()
static ::ExceptionBase & ExcNotImplemented()
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
Definition: polynomial.cc:823
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
Definition: polynomial.cc:1230