Reference documentation for deal.II version Git c064f80013 2020-07-04 08:24:28 -0400
\(\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 
100  number
101  value(const number x) const;
102 
113  void
114  value(const number x, std::vector<number> &values) const;
115 
127  void
128  value(const number x,
129  const unsigned int n_derivatives,
130  number * values) const;
131 
137  unsigned int
138  degree() const;
139 
147  void
148  scale(const number factor);
149 
165  template <typename number2>
166  void
167  shift(const number2 offset);
168 
173  derivative() const;
174 
180  primitive() const;
181 
186  operator*=(const double s);
187 
192  operator*=(const Polynomial<number> &p);
193 
198  operator+=(const Polynomial<number> &p);
199 
204  operator-=(const Polynomial<number> &p);
205 
209  bool
210  operator==(const Polynomial<number> &p) const;
211 
215  void
216  print(std::ostream &out) const;
217 
222  template <class Archive>
223  void
224  serialize(Archive &ar, const unsigned int version);
225 
229  virtual std::size_t
230  memory_consumption() const;
231 
232  protected:
236  static void
237  scale(std::vector<number> &coefficients, const number factor);
238 
242  template <typename number2>
243  static void
244  shift(std::vector<number> &coefficients, const number2 shift);
245 
249  static void
250  multiply(std::vector<number> &coefficients, const number factor);
251 
257  void
259 
268  std::vector<number> coefficients;
269 
275 
280  std::vector<number> lagrange_support_points;
281 
287  };
288 
289 
294  template <typename number>
295  class Monomial : public Polynomial<number>
296  {
297  public:
302  Monomial(const unsigned int n, const double coefficient = 1.);
303 
310  static std::vector<Polynomial<number>>
311  generate_complete_basis(const unsigned int degree);
312 
313  private:
317  static std::vector<number>
318  make_vector(unsigned int n, const double coefficient);
319  };
320 
321 
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 
395  class Legendre : public Polynomial<double>
396  {
397  public:
401  Legendre(const unsigned int p);
402 
409  static std::vector<Polynomial<double>>
410  generate_complete_basis(const unsigned int degree);
411  };
412 
432  class Lobatto : public Polynomial<double>
433  {
434  public:
439  Lobatto(const unsigned int p = 0);
440 
445  static std::vector<Polynomial<double>>
446  generate_complete_basis(const unsigned int p);
447 
448  private:
452  std::vector<double>
453  compute_coefficients(const unsigned int p);
454  };
455 
456 
457 
495  class Hierarchical : public Polynomial<double>
496  {
497  public:
502  Hierarchical(const unsigned int p);
503 
514  static std::vector<Polynomial<double>>
515  generate_complete_basis(const unsigned int degree);
516 
517  private:
521  static void
522  compute_coefficients(const unsigned int p);
523 
528  static const std::vector<double> &
529  get_coefficients(const unsigned int p);
530 
539  static std::vector<std::unique_ptr<const std::vector<double>>>
541  };
542 
543 
544 
572  class HermiteInterpolation : public Polynomial<double>
573  {
574  public:
579  HermiteInterpolation(const unsigned int p);
580 
586  static std::vector<Polynomial<double>>
587  generate_complete_basis(const unsigned int p);
588  };
589 
590 
591 
693  class HermiteLikeInterpolation : public Polynomial<double>
694  {
695  public:
700  HermiteLikeInterpolation(const unsigned int degree,
701  const unsigned int index);
702 
707  static std::vector<Polynomial<double>>
708  generate_complete_basis(const unsigned int degree);
709  };
710 
711 
712 
713  /*
714  * Evaluate a Jacobi polynomial @f$ P_n^{\alpha, \beta}(x) @f$ specified by the
715  * parameters @p alpha, @p beta, @p n, where @p n is the degree of the
716  * Jacobi polynomial.
717  *
718  * @note The Jacobi polynomials are not orthonormal and are defined on the
719  * unit interval @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often
720  * used in literature. @p x is the point of evaluation.
721  */
722  template <typename Number>
723  Number
724  jacobi_polynomial_value(const unsigned int degree,
725  const int alpha,
726  const int beta,
727  const Number x);
728 
729 
742  template <typename Number>
743  std::vector<Number>
744  jacobi_polynomial_roots(const unsigned int degree,
745  const int alpha,
746  const int beta);
747 } // namespace Polynomials
748 
749 
752 /* -------------------------- inline functions --------------------- */
753 
754 namespace Polynomials
755 {
756  template <typename number>
758  : in_lagrange_product_form(false)
759  , lagrange_weight(1.)
760  {}
761 
762 
763 
764  template <typename number>
765  inline unsigned int
767  {
768  if (in_lagrange_product_form == true)
769  {
770  return lagrange_support_points.size();
771  }
772  else
773  {
774  Assert(coefficients.size() > 0, ExcEmptyObject());
775  return coefficients.size() - 1;
776  }
777  }
778 
779 
780 
781  template <typename number>
782  inline number
783  Polynomial<number>::value(const number x) const
784  {
785  if (in_lagrange_product_form == false)
786  {
787  Assert(coefficients.size() > 0, ExcEmptyObject());
788 
789  // Horner scheme
790  const unsigned int m = coefficients.size();
791  number value = coefficients.back();
792  for (int k = m - 2; k >= 0; --k)
793  value = value * x + coefficients[k];
794  return value;
795  }
796  else
797  {
798  // direct evaluation of Lagrange polynomial
799  const unsigned int m = lagrange_support_points.size();
800  number value = 1.;
801  for (unsigned int j = 0; j < m; ++j)
802  value *= x - lagrange_support_points[j];
803  value *= lagrange_weight;
804  return value;
805  }
806  }
807 
808 
809 
810  template <typename number>
811  template <class Archive>
812  inline void
813  Polynomial<number>::serialize(Archive &ar, const unsigned int)
814  {
815  // forward to serialization function in the base class.
816  ar &static_cast<Subscriptor &>(*this);
817  ar &coefficients;
818  ar &in_lagrange_product_form;
819  ar &lagrange_support_points;
820  ar &lagrange_weight;
821  }
822 
823 
824 
825  template <typename Number>
826  Number
827  jacobi_polynomial_value(const unsigned int degree,
828  const int alpha,
829  const int beta,
830  const Number x)
831  {
832  Assert(alpha >= 0 && beta >= 0,
833  ExcNotImplemented("Negative alpha/beta coefficients not supported"));
834  // the Jacobi polynomial is evaluated using a recursion formula.
835  Number p0, p1;
836 
837  // The recursion formula is defined for the interval [-1, 1], so rescale
838  // to that interval here
839  const Number xeval = Number(-1) + 2. * x;
840 
841  // initial values P_0(x), P_1(x):
842  p0 = 1.0;
843  if (degree == 0)
844  return p0;
845  p1 = ((alpha + beta + 2) * xeval + (alpha - beta)) / 2;
846  if (degree == 1)
847  return p1;
848 
849  for (unsigned int i = 1; i < degree; ++i)
850  {
851  const Number v = 2 * i + (alpha + beta);
852  const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
853  const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
854  const Number a3 = v * (v + 1) * (v + 2);
855  const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
856 
857  const Number pn = ((a2 + a3 * xeval) * p1 - a4 * p0) / a1;
858  p0 = p1;
859  p1 = pn;
860  }
861  return p1;
862  }
863 
864 
865 
866  template <typename Number>
867  std::vector<Number>
868  jacobi_polynomial_roots(const unsigned int degree,
869  const int alpha,
870  const int beta)
871  {
872  std::vector<Number> x(degree, 0.5);
873 
874  // compute zeros with a Newton algorithm.
875 
876  // Set tolerance. For long double we might not always get the additional
877  // precision in a run time environment (e.g. with valgrind), so we must
878  // limit the tolerance to double. Since we do a Newton iteration, doing
879  // one more iteration after the residual has indicated convergence will be
880  // enough for all number types due to the quadratic convergence of
881  // Newton's method
882 
883  const Number tolerance =
884  4 * std::max(static_cast<Number>(std::numeric_limits<double>::epsilon()),
886 
887  // The following implementation follows closely the one given in the
888  // appendix of the book by Karniadakis and Sherwin: Spectral/hp element
889  // methods for computational fluid dynamics (Oxford University Press,
890  // 2005)
891 
892  // If symmetric, we only need to compute the half of points
893  const unsigned int n_points = (alpha == beta ? degree / 2 : degree);
894  for (unsigned int k = 0; k < n_points; ++k)
895  {
896  // we take the zeros of the Chebyshev polynomial (alpha=beta=-0.5) as
897  // initial values, corrected by the initial value
898  Number r = 0.5 - 0.5 * std::cos(static_cast<Number>(2 * k + 1) /
899  (2 * degree) * numbers::PI);
900  if (k > 0)
901  r = (r + x[k - 1]) / 2;
902 
903  unsigned int converged = numbers::invalid_unsigned_int;
904  for (unsigned int it = 1; it < 1000; ++it)
905  {
906  Number s = 0.;
907  for (unsigned int i = 0; i < k; ++i)
908  s += 1. / (r - x[i]);
909 
910  // derivative of P_n^{alpha,beta}, rescaled to [0, 1]
911  const Number J_x =
912  (alpha + beta + degree + 1) *
913  jacobi_polynomial_value(degree - 1, alpha + 1, beta + 1, r);
914 
915  // value of P_n^{alpha,beta}
916  const Number f = jacobi_polynomial_value(degree, alpha, beta, r);
917  const Number delta = f / (f * s - J_x);
918  r += delta;
919  if (converged == numbers::invalid_unsigned_int &&
920  std::abs(delta) < tolerance)
921  converged = it;
922 
923  // do one more iteration to ensure accuracy also for tighter
924  // types than double (e.g. long double)
925  if (it == converged + 1)
926  break;
927  }
928 
930  ExcMessage("Newton iteration for zero of Jacobi polynomial "
931  "did not converge."));
932 
933  x[k] = r;
934  }
935 
936  // in case we assumed symmetry, fill up the missing values
937  for (unsigned int k = n_points; k < degree; ++k)
938  x[k] = 1.0 - x[degree - k - 1];
939 
940  return x;
941  }
942 
943 } // namespace Polynomials
945 
946 #endif
static const unsigned int invalid_unsigned_int
Definition: types.h:191
void serialize(Archive &ar, const unsigned int version)
Definition: polynomial.h:813
Polynomial< number > & operator*=(const double s)
Definition: polynomial.cc:336
void scale(const number factor)
Definition: polynomial.cc:298
void transform_into_standard_form()
Definition: polynomial.cc:244
bool operator==(const Polynomial< number > &p) const
Definition: polynomial.cc:479
Polynomial< number > derivative() const
Definition: polynomial.cc:591
void shift(const number2 offset)
Definition: polynomial.cc:572
unsigned int degree() const
Definition: polynomial.h:766
Polynomial< number > & operator-=(const Polynomial< number > &p)
Definition: polynomial.cc:443
static std::vector< std::unique_ptr< const std::vector< double > > > recursive_coefficients
Definition: polynomial.h:540
std::vector< Number > jacobi_polynomial_roots(const unsigned int degree, const int alpha, const int beta)
Definition: polynomial.h:868
number value(const number x) const
Definition: polynomial.h:783
Definition: point.h:110
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
void print(std::ostream &out) const
Definition: polynomial.cc:647
Number jacobi_polynomial_value(const unsigned int degree, const int alpha, const int beta, const Number x)
Definition: polynomial.h:827
static ::ExceptionBase & ExcMessage(std::string arg1)
static void multiply(std::vector< number > &coefficients, const number factor)
Definition: polynomial.cc:323
#define Assert(cond, exc)
Definition: exceptions.h:1403
Polynomial< number > & operator+=(const Polynomial< number > &p)
Definition: polynomial.cc:401
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
virtual std::size_t memory_consumption() const
Definition: polynomial.cc:666
std::vector< number > coefficients
Definition: polynomial.h:268
Polynomial< number > primitive() const
Definition: polynomial.cc:620
std::vector< number > lagrange_support_points
Definition: polynomial.h:280
static constexpr double PI
Definition: numbers.h:231
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
static ::ExceptionBase & ExcEmptyObject()
static ::ExceptionBase & ExcNotImplemented()
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:834