Reference documentation for deal.II version Git 3ff8e288d2 2020-02-20 17:30:00 +0100
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
231  virtual std::size_t
232  memory_consumption() const;
233 
234  protected:
238  static void
239  scale(std::vector<number> &coefficients, const number factor);
240 
244  template <typename number2>
245  static void
246  shift(std::vector<number> &coefficients, const number2 shift);
247 
251  static void
252  multiply(std::vector<number> &coefficients, const number factor);
253 
259  void
261 
270  std::vector<number> coefficients;
271 
277 
282  std::vector<number> lagrange_support_points;
283 
289  };
290 
291 
298  template <typename number>
299  class Monomial : public Polynomial<number>
300  {
301  public:
306  Monomial(const unsigned int n, const double coefficient = 1.);
307 
314  static std::vector<Polynomial<number>>
315  generate_complete_basis(const unsigned int degree);
316 
317  private:
321  static std::vector<number>
322  make_vector(unsigned int n, const double coefficient);
323  };
324 
325 
343  class LagrangeEquidistant : public Polynomial<double>
344  {
345  public:
351  LagrangeEquidistant(const unsigned int n, const unsigned int support_point);
352 
361  static std::vector<Polynomial<double>>
362  generate_complete_basis(const unsigned int degree);
363 
364  private:
369  static void
370  compute_coefficients(const unsigned int n,
371  const unsigned int support_point,
372  std::vector<double> &a);
373  };
374 
375 
376 
383  std::vector<Polynomial<double>>
384  generate_complete_Lagrange_basis(const std::vector<Point<1>> &points);
385 
386 
387 
403  class Legendre : public Polynomial<double>
404  {
405  public:
409  Legendre(const unsigned int p);
410 
417  static std::vector<Polynomial<double>>
418  generate_complete_basis(const unsigned int degree);
419  };
420 
442  class Lobatto : public Polynomial<double>
443  {
444  public:
449  Lobatto(const unsigned int p = 0);
450 
455  static std::vector<Polynomial<double>>
456  generate_complete_basis(const unsigned int p);
457 
458  private:
462  std::vector<double>
463  compute_coefficients(const unsigned int p);
464  };
465 
466 
467 
507  class Hierarchical : public Polynomial<double>
508  {
509  public:
514  Hierarchical(const unsigned int p);
515 
526  static std::vector<Polynomial<double>>
527  generate_complete_basis(const unsigned int degree);
528 
529  private:
533  static void
534  compute_coefficients(const unsigned int p);
535 
540  static const std::vector<double> &
541  get_coefficients(const unsigned int p);
542 
551  static std::vector<std::unique_ptr<const std::vector<double>>>
553  };
554 
555 
556 
587  class HermiteInterpolation : public Polynomial<double>
588  {
589  public:
594  HermiteInterpolation(const unsigned int p);
595 
601  static std::vector<Polynomial<double>>
602  generate_complete_basis(const unsigned int p);
603  };
604 
605 
606 
711  class HermiteLikeInterpolation : public Polynomial<double>
712  {
713  public:
718  HermiteLikeInterpolation(const unsigned int degree,
719  const unsigned int index);
720 
725  static std::vector<Polynomial<double>>
726  generate_complete_basis(const unsigned int degree);
727  };
728 
729 
730 
731  /*
732  * Evaluate a Jacobi polynomial @f$ P_n^{\alpha, \beta}(x) @f$ specified by the
733  * parameters @p alpha, @p beta, @p n, where @p n is the degree of the
734  * Jacobi polynomial.
735  *
736  * @note The Jacobi polynomials are not orthonormal and are defined on the
737  * unit interval @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often
738  * used in literature. @p x is the point of evaluation.
739  */
740  template <typename Number>
741  Number
742  jacobi_polynomial_value(const unsigned int degree,
743  const int alpha,
744  const int beta,
745  const Number x);
746 
747 
760  template <typename Number>
761  std::vector<Number>
762  jacobi_polynomial_roots(const unsigned int degree,
763  const int alpha,
764  const int beta);
765 } // namespace Polynomials
766 
767 
770 /* -------------------------- inline functions --------------------- */
771 
772 namespace Polynomials
773 {
774  template <typename number>
776  : in_lagrange_product_form(false)
777  , lagrange_weight(1.)
778  {}
779 
780 
781 
782  template <typename number>
783  inline unsigned int
785  {
786  if (in_lagrange_product_form == true)
787  {
788  return lagrange_support_points.size();
789  }
790  else
791  {
792  Assert(coefficients.size() > 0, ExcEmptyObject());
793  return coefficients.size() - 1;
794  }
795  }
796 
797 
798 
799  template <typename number>
800  inline number
801  Polynomial<number>::value(const number x) const
802  {
803  if (in_lagrange_product_form == false)
804  {
805  Assert(coefficients.size() > 0, ExcEmptyObject());
806 
807  // Horner scheme
808  const unsigned int m = coefficients.size();
809  number value = coefficients.back();
810  for (int k = m - 2; k >= 0; --k)
811  value = value * x + coefficients[k];
812  return value;
813  }
814  else
815  {
816  // direct evaluation of Lagrange polynomial
817  const unsigned int m = lagrange_support_points.size();
818  number value = 1.;
819  for (unsigned int j = 0; j < m; ++j)
820  value *= x - lagrange_support_points[j];
821  value *= lagrange_weight;
822  return value;
823  }
824  }
825 
826 
827 
828  template <typename number>
829  template <class Archive>
830  inline void
831  Polynomial<number>::serialize(Archive &ar, const unsigned int)
832  {
833  // forward to serialization function in the base class.
834  ar &static_cast<Subscriptor &>(*this);
835  ar &coefficients;
836  ar &in_lagrange_product_form;
837  ar &lagrange_support_points;
838  ar &lagrange_weight;
839  }
840 
841 
842 
843  template <typename Number>
844  Number
845  jacobi_polynomial_value(const unsigned int degree,
846  const int alpha,
847  const int beta,
848  const Number x)
849  {
850  Assert(alpha >= 0 && beta >= 0,
851  ExcNotImplemented("Negative alpha/beta coefficients not supported"));
852  // the Jacobi polynomial is evaluated using a recursion formula.
853  Number p0, p1;
854 
855  // The recursion formula is defined for the interval [-1, 1], so rescale
856  // to that interval here
857  const Number xeval = Number(-1) + 2. * x;
858 
859  // initial values P_0(x), P_1(x):
860  p0 = 1.0;
861  if (degree == 0)
862  return p0;
863  p1 = ((alpha + beta + 2) * xeval + (alpha - beta)) / 2;
864  if (degree == 1)
865  return p1;
866 
867  for (unsigned int i = 1; i < degree; ++i)
868  {
869  const Number v = 2 * i + (alpha + beta);
870  const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
871  const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
872  const Number a3 = v * (v + 1) * (v + 2);
873  const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
874 
875  const Number pn = ((a2 + a3 * xeval) * p1 - a4 * p0) / a1;
876  p0 = p1;
877  p1 = pn;
878  }
879  return p1;
880  }
881 
882 
883 
884  template <typename Number>
885  std::vector<Number>
886  jacobi_polynomial_roots(const unsigned int degree,
887  const int alpha,
888  const int beta)
889  {
890  std::vector<Number> x(degree, 0.5);
891 
892  // compute zeros with a Newton algorithm.
893 
894  // Set tolerance. For long double we might not always get the additional
895  // precision in a run time environment (e.g. with valgrind), so we must
896  // limit the tolerance to double. Since we do a Newton iteration, doing
897  // one more iteration after the residual has indicated convergence will be
898  // enough for all number types due to the quadratic convergence of
899  // Newton's method
900 
901  const Number tolerance =
902  4 * std::max(static_cast<Number>(std::numeric_limits<double>::epsilon()),
903  std::numeric_limits<Number>::epsilon());
904 
905  // The following implementation follows closely the one given in the
906  // appendix of the book by Karniadakis and Sherwin: Spectral/hp element
907  // methods for computational fluid dynamics (Oxford University Press,
908  // 2005)
909 
910  // If symmetric, we only need to compute the half of points
911  const unsigned int n_points = (alpha == beta ? degree / 2 : degree);
912  for (unsigned int k = 0; k < n_points; ++k)
913  {
914  // we take the zeros of the Chebyshev polynomial (alpha=beta=-0.5) as
915  // initial values, corrected by the initial value
916  Number r = 0.5 - 0.5 * std::cos(static_cast<Number>(2 * k + 1) /
917  (2 * degree) * numbers::PI);
918  if (k > 0)
919  r = (r + x[k - 1]) / 2;
920 
921  unsigned int converged = numbers::invalid_unsigned_int;
922  for (unsigned int it = 1; it < 1000; ++it)
923  {
924  Number s = 0.;
925  for (unsigned int i = 0; i < k; ++i)
926  s += 1. / (r - x[i]);
927 
928  // derivative of P_n^{alpha,beta}, rescaled to [0, 1]
929  const Number J_x =
930  (alpha + beta + degree + 1) *
931  jacobi_polynomial_value(degree - 1, alpha + 1, beta + 1, r);
932 
933  // value of P_n^{alpha,beta}
934  const Number f = jacobi_polynomial_value(degree, alpha, beta, r);
935  const Number delta = f / (f * s - J_x);
936  r += delta;
937  if (converged == numbers::invalid_unsigned_int &&
938  std::abs(delta) < tolerance)
939  converged = it;
940 
941  // do one more iteration to ensure accuracy also for tighter
942  // types than double (e.g. long double)
943  if (it == converged + 1)
944  break;
945  }
946 
948  ExcMessage("Newton iteration for zero of Jacobi polynomial "
949  "did not converge."));
950 
951  x[k] = r;
952  }
953 
954  // in case we assumed symmetry, fill up the missing values
955  for (unsigned int k = n_points; k < degree; ++k)
956  x[k] = 1.0 - x[degree - k - 1];
957 
958  return x;
959  }
960 
961 } // namespace Polynomials
962 DEAL_II_NAMESPACE_CLOSE
963 
964 #endif
static const unsigned int invalid_unsigned_int
Definition: types.h:187
void serialize(Archive &ar, const unsigned int version)
Definition: polynomial.h:831
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:784
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:552
std::vector< Number > jacobi_polynomial_roots(const unsigned int degree, const int alpha, const int beta)
Definition: polynomial.h:886
number value(const number x) const
Definition: polynomial.h:801
Definition: point.h:110
void print(std::ostream &out) const
Definition: polynomial.cc:647
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:1411
Polynomial< number > & operator+=(const Polynomial< number > &p)
Definition: polynomial.cc:401
virtual std::size_t memory_consumption() const
Definition: polynomial.cc:666
std::vector< number > coefficients
Definition: polynomial.h:270
Polynomial< number > primitive() const
Definition: polynomial.cc:620
std::vector< number > lagrange_support_points
Definition: polynomial.h:282
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:834