Reference documentation for deal.II version Git 7e76c9f1da 2019-09-19 09:07:25 +0200
\(\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 
705  class HermiteLikeInterpolation : public Polynomial<double>
706  {
707  public:
712  HermiteLikeInterpolation(const unsigned int degree,
713  const unsigned int index);
714 
719  static std::vector<Polynomial<double>>
720  generate_complete_basis(const unsigned int degree);
721  };
722 
723 
724 
725  /*
726  * Evaluate a Jacobi polynomial @f$ P_n^{\alpha, \beta}(x) @f$ specified by the
727  * parameters @p alpha, @p beta, @p n, where @p n is the degree of the
728  * Jacobi polynomial.
729  *
730  * @note The Jacobi polynomials are not orthonormal and are defined on the
731  * unit interval @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often
732  * used in literature. @p x is the point of evaluation.
733  */
734  template <typename Number>
735  Number
736  jacobi_polynomial_value(const unsigned int degree,
737  const int alpha,
738  const int beta,
739  const Number x);
740 
741 
754  template <typename Number>
755  std::vector<Number>
756  jacobi_polynomial_roots(const unsigned int degree,
757  const int alpha,
758  const int beta);
759 } // namespace Polynomials
760 
761 
764 /* -------------------------- inline functions --------------------- */
765 
766 namespace Polynomials
767 {
768  template <typename number>
770  : in_lagrange_product_form(false)
771  , lagrange_weight(1.)
772  {}
773 
774 
775 
776  template <typename number>
777  inline unsigned int
779  {
780  if (in_lagrange_product_form == true)
781  {
782  return lagrange_support_points.size();
783  }
784  else
785  {
786  Assert(coefficients.size() > 0, ExcEmptyObject());
787  return coefficients.size() - 1;
788  }
789  }
790 
791 
792 
793  template <typename number>
794  inline number
795  Polynomial<number>::value(const number x) const
796  {
797  if (in_lagrange_product_form == false)
798  {
799  Assert(coefficients.size() > 0, ExcEmptyObject());
800 
801  // Horner scheme
802  const unsigned int m = coefficients.size();
803  number value = coefficients.back();
804  for (int k = m - 2; k >= 0; --k)
805  value = value * x + coefficients[k];
806  return value;
807  }
808  else
809  {
810  // direct evaluation of Lagrange polynomial
811  const unsigned int m = lagrange_support_points.size();
812  number value = 1.;
813  for (unsigned int j = 0; j < m; ++j)
814  value *= x - lagrange_support_points[j];
815  value *= lagrange_weight;
816  return value;
817  }
818  }
819 
820 
821 
822  template <typename number>
823  template <class Archive>
824  inline void
825  Polynomial<number>::serialize(Archive &ar, const unsigned int)
826  {
827  // forward to serialization function in the base class.
828  ar &static_cast<Subscriptor &>(*this);
829  ar &coefficients;
830  ar &in_lagrange_product_form;
831  ar &lagrange_support_points;
832  ar &lagrange_weight;
833  }
834 
835 
836 
837  template <typename Number>
838  Number
839  jacobi_polynomial_value(const unsigned int degree,
840  const int alpha,
841  const int beta,
842  const Number x)
843  {
844  Assert(alpha >= 0 && beta >= 0,
845  ExcNotImplemented("Negative alpha/beta coefficients not supported"));
846  // the Jacobi polynomial is evaluated using a recursion formula.
847  Number p0, p1;
848 
849  // The recursion formula is defined for the interval [-1, 1], so rescale
850  // to that interval here
851  const Number xeval = Number(-1) + 2. * x;
852 
853  // initial values P_0(x), P_1(x):
854  p0 = 1.0;
855  if (degree == 0)
856  return p0;
857  p1 = ((alpha + beta + 2) * xeval + (alpha - beta)) / 2;
858  if (degree == 1)
859  return p1;
860 
861  for (unsigned int i = 1; i < degree; ++i)
862  {
863  const Number v = 2 * i + (alpha + beta);
864  const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
865  const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
866  const Number a3 = v * (v + 1) * (v + 2);
867  const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
868 
869  const Number pn = ((a2 + a3 * xeval) * p1 - a4 * p0) / a1;
870  p0 = p1;
871  p1 = pn;
872  }
873  return p1;
874  }
875 
876 
877 
878  template <typename Number>
879  std::vector<Number>
880  jacobi_polynomial_roots(const unsigned int degree,
881  const int alpha,
882  const int beta)
883  {
884  std::vector<Number> x(degree, 0.5);
885 
886  // compute zeros with a Newton algorithm.
887 
888  // Set tolerance. For long double we might not always get the additional
889  // precision in a run time environment (e.g. with valgrind), so we must
890  // limit the tolerance to double. Since we do a Newton iteration, doing
891  // one more iteration after the residual has indicated convergence will be
892  // enough for all number types due to the quadratic convergence of
893  // Newton's method
894 
895  const Number tolerance =
896  4 * std::max(static_cast<Number>(std::numeric_limits<double>::epsilon()),
897  std::numeric_limits<Number>::epsilon());
898 
899  // The following implementation follows closely the one given in the
900  // appendix of the book by Karniadakis and Sherwin: Spectral/hp element
901  // methods for computational fluid dynamics (Oxford University Press,
902  // 2005)
903 
904  // If symmetric, we only need to compute the half of points
905  const unsigned int n_points = (alpha == beta ? degree / 2 : degree);
906  for (unsigned int k = 0; k < n_points; ++k)
907  {
908  // we take the zeros of the Chebyshev polynomial (alpha=beta=-0.5) as
909  // initial values, corrected by the initial value
910  Number r = 0.5 - 0.5 * std::cos(static_cast<Number>(2 * k + 1) /
911  (2 * degree) * numbers::PI);
912  if (k > 0)
913  r = (r + x[k - 1]) / 2;
914 
915  unsigned int converged = numbers::invalid_unsigned_int;
916  for (unsigned int it = 1; it < 1000; ++it)
917  {
918  Number s = 0.;
919  for (unsigned int i = 0; i < k; ++i)
920  s += 1. / (r - x[i]);
921 
922  // derivative of P_n^{alpha,beta}, rescaled to [0, 1]
923  const Number J_x =
924  (alpha + beta + degree + 1) *
925  jacobi_polynomial_value(degree - 1, alpha + 1, beta + 1, r);
926 
927  // value of P_n^{alpha,beta}
928  const Number f = jacobi_polynomial_value(degree, alpha, beta, r);
929  const Number delta = f / (f * s - J_x);
930  r += delta;
931  if (converged == numbers::invalid_unsigned_int &&
932  std::abs(delta) < tolerance)
933  converged = it;
934 
935  // do one more iteration to ensure accuracy also for tighter
936  // types than double (e.g. long double)
937  if (it == converged + 1)
938  break;
939  }
940 
942  ExcMessage("Newton iteration for zero of Jacobi polynomial "
943  "did not converge."));
944 
945  x[k] = r;
946  }
947 
948  // in case we assumed symmetry, fill up the missing values
949  for (unsigned int k = n_points; k < degree; ++k)
950  x[k] = 1.0 - x[degree - k - 1];
951 
952  return x;
953  }
954 
955 } // namespace Polynomials
956 DEAL_II_NAMESPACE_CLOSE
957 
958 #endif
static const unsigned int invalid_unsigned_int
Definition: types.h:187
void serialize(Archive &ar, const unsigned int version)
Definition: polynomial.h:825
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
Polynomial< number > derivative() const
Definition: polynomial.cc:590
void shift(const number2 offset)
Definition: polynomial.cc:571
unsigned int degree() const
Definition: polynomial.h:778
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:880
number value(const number x) const
Definition: polynomial.h:795
Definition: point.h:110
void print(std::ostream &out) const
Definition: polynomial.cc:646
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
std::vector< number > coefficients
Definition: polynomial.h:264
Polynomial< number > primitive() const
Definition: polynomial.cc:619
std::vector< number > lagrange_support_points
Definition: polynomial.h:276
static constexpr double PI
Definition: numbers.h:238
static ::ExceptionBase & ExcEmptyObject()
static ::ExceptionBase & ExcNotImplemented()
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
Definition: polynomial.cc:823