Reference documentation for deal.II version GIT relicensing-43-g3fc2fae7c3 2024-03-03 21:10:02+00:00
\(\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\}}\)
Loading...
Searching...
No Matches
polynomial.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2000 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_polynomial_h
16#define dealii_polynomial_h
17
18
19
20#include <deal.II/base/config.h>
21
23#include <deal.II/base/point.h>
25
26#include <array>
27#include <limits>
28#include <memory>
29#include <shared_mutex>
30#include <vector>
31
33
43namespace Polynomials
44{
64 template <typename number>
65 class Polynomial : public Subscriptor
66 {
67 public:
76 Polynomial(const std::vector<number> &coefficients);
77
81 Polynomial(const unsigned int n);
82
91 const unsigned int evaluation_point);
92
97
107 number
108 value(const number x) const;
109
120 void
121 value(const number x, std::vector<number> &values) const;
122
141 template <typename Number2>
142 void
143 value(const Number2 x,
144 const unsigned int n_derivatives,
145 Number2 *values) const;
146
159 template <std::size_t n_entries, typename Number2>
160 void
161 values_of_array(const std::array<Number2, n_entries> &points,
162 const unsigned int n_derivatives,
163 std::array<Number2, n_entries> *values) const;
164
170 unsigned int
171 degree() const;
172
180 void
181 scale(const number factor);
182
198 template <typename number2>
199 void
200 shift(const number2 offset);
201
206 derivative() const;
207
213 primitive() const;
214
219 operator*=(const double s);
220
226
232
238
242 bool
243 operator==(const Polynomial<number> &p) const;
244
248 void
249 print(std::ostream &out) const;
250
256 template <class Archive>
257 void
258 serialize(Archive &ar, const unsigned int version);
259
263 virtual std::size_t
264 memory_consumption() const;
265
266 protected:
270 static void
271 scale(std::vector<number> &coefficients, const number factor);
272
276 template <typename number2>
277 static void
278 shift(std::vector<number> &coefficients, const number2 shift);
279
283 static void
284 multiply(std::vector<number> &coefficients, const number factor);
285
291 void
293
302 std::vector<number> coefficients;
303
309
314 std::vector<number> lagrange_support_points;
315
321 };
322
323
328 template <typename number>
329 class Monomial : public Polynomial<number>
330 {
331 public:
336 Monomial(const unsigned int n, const double coefficient = 1.);
337
344 static std::vector<Polynomial<number>>
345 generate_complete_basis(const unsigned int degree);
346
347 private:
351 static std::vector<number>
352 make_vector(unsigned int n, const double coefficient);
353 };
354
355
371 class LagrangeEquidistant : public Polynomial<double>
372 {
373 public:
379 LagrangeEquidistant(const unsigned int n, const unsigned int support_point);
380
389 static std::vector<Polynomial<double>>
390 generate_complete_basis(const unsigned int degree);
391
392 private:
397 static void
398 compute_coefficients(const unsigned int n,
399 const unsigned int support_point,
400 std::vector<double> &a);
401 };
402
403
404
411 std::vector<Polynomial<double>>
412 generate_complete_Lagrange_basis(const std::vector<Point<1>> &points);
413
414
415
429 class Legendre : public Polynomial<double>
430 {
431 public:
435 Legendre(const unsigned int p);
436
443 static std::vector<Polynomial<double>>
444 generate_complete_basis(const unsigned int degree);
445 };
446
466 class Lobatto : public Polynomial<double>
467 {
468 public:
473 Lobatto(const unsigned int p = 0);
474
479 static std::vector<Polynomial<double>>
480 generate_complete_basis(const unsigned int p);
481
482 private:
486 std::vector<double>
487 compute_coefficients(const unsigned int p);
488 };
489
490
491
529 class Hierarchical : public Polynomial<double>
530 {
531 public:
536 Hierarchical(const unsigned int p);
537
548 static std::vector<Polynomial<double>>
549 generate_complete_basis(const unsigned int degree);
550
551 private:
555 static void
556 compute_coefficients(const unsigned int p);
557
562 static const std::vector<double> &
563 get_coefficients(const unsigned int p);
564
573 static std::vector<std::unique_ptr<const std::vector<double>>>
575
580 static std::shared_mutex coefficients_lock;
581 };
582
583
584
612 class HermiteInterpolation : public Polynomial<double>
613 {
614 public:
619 HermiteInterpolation(const unsigned int p);
620
626 static std::vector<Polynomial<double>>
627 generate_complete_basis(const unsigned int p);
628 };
629
630
631
734 {
735 public:
740 HermiteLikeInterpolation(const unsigned int degree,
741 const unsigned int index);
742
747 static std::vector<Polynomial<double>>
748 generate_complete_basis(const unsigned int degree);
749 };
750
751
752
753 /*
754 * Evaluate a Jacobi polynomial @f$ P_n^{\alpha, \beta}(x) @f$ specified by the
755 * parameters @p alpha, @p beta, @p n, where @p n is the degree of the
756 * Jacobi polynomial.
757 *
758 * @note The Jacobi polynomials are not orthonormal and are defined on the
759 * unit interval @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often
760 * used in literature. @p x is the point of evaluation.
761 */
762 template <typename Number>
763 Number
764 jacobi_polynomial_value(const unsigned int degree,
765 const int alpha,
766 const int beta,
767 const Number x);
768
769
782 template <typename Number>
783 std::vector<Number>
784 jacobi_polynomial_roots(const unsigned int degree,
785 const int alpha,
786 const int beta);
787} // namespace Polynomials
788
789
792/* -------------------------- inline functions --------------------- */
793
794namespace Polynomials
795{
796 template <typename number>
798 : in_lagrange_product_form(false)
799 , lagrange_weight(1.)
800 {}
801
802
803
804 template <typename number>
805 inline unsigned int
807 {
808 if (in_lagrange_product_form == true)
809 {
810 return lagrange_support_points.size();
811 }
812 else
813 {
814 Assert(coefficients.size() > 0, ExcEmptyObject());
815 return coefficients.size() - 1;
816 }
817 }
818
819
820
821 template <typename number>
822 inline number
823 Polynomial<number>::value(const number x) const
824 {
825 if (in_lagrange_product_form == false)
826 {
827 Assert(coefficients.size() > 0, ExcEmptyObject());
828
829 // Horner scheme
830 const unsigned int m = coefficients.size();
831 number value = coefficients.back();
832 for (int k = m - 2; k >= 0; --k)
833 value = value * x + coefficients[k];
834 return value;
835 }
836 else
837 {
838 // direct evaluation of Lagrange polynomial
839 const unsigned int m = lagrange_support_points.size();
840 number value = 1.;
841 for (unsigned int j = 0; j < m; ++j)
842 value *= x - lagrange_support_points[j];
843 value *= lagrange_weight;
844 return value;
845 }
846 }
847
848
849
850 template <typename number>
851 template <typename Number2>
852 inline void
854 const unsigned int n_derivatives,
855 Number2 *values) const
856 {
857 values_of_array(std::array<Number2, 1ul>{{x}},
858 n_derivatives,
859 reinterpret_cast<std::array<Number2, 1ul> *>(values));
860 }
861
862
863
864 template <typename number>
865 template <std::size_t n_entries, typename Number2>
866 inline void
868 const std::array<Number2, n_entries> &x,
869 const unsigned int n_derivatives,
870 std::array<Number2, n_entries> *values) const
871 {
872 // evaluate Lagrange polynomial and derivatives
873 if (in_lagrange_product_form == true)
874 {
875 // to compute the value and all derivatives of a polynomial of the
876 // form (x-x_1)*(x-x_2)*...*(x-x_n), expand the derivatives like
877 // automatic differentiation does.
878 const unsigned int n_supp = lagrange_support_points.size();
879 const number weight = lagrange_weight;
880 switch (n_derivatives)
881 {
882 default:
883 for (unsigned int e = 0; e < n_entries; ++e)
884 values[0][e] = weight;
885 for (unsigned int k = 1; k <= n_derivatives; ++k)
886 for (unsigned int e = 0; e < n_entries; ++e)
887 values[k][e] = 0.;
888 for (unsigned int i = 0; i < n_supp; ++i)
889 {
890 std::array<Number2, n_entries> v = x;
891 for (unsigned int e = 0; e < n_entries; ++e)
892 v[e] -= lagrange_support_points[i];
893
894 // multiply by (x-x_i) and compute action on all derivatives,
895 // too (inspired from automatic differentiation: implement the
896 // product rule for the old value and the new variable 'v',
897 // i.e., expand value v and derivative one). since we reuse a
898 // value from the next lower derivative from the steps before,
899 // need to start from the highest derivative
900 for (unsigned int k = n_derivatives; k > 0; --k)
901 for (unsigned int e = 0; e < n_entries; ++e)
902 values[k][e] = (values[k][e] * v[e] + values[k - 1][e]);
903 for (unsigned int e = 0; e < n_entries; ++e)
904 values[0][e] *= v[e];
905 }
906 // finally, multiply derivatives by k! to transform the product
907 // p_n = p^(n)(x)/k! into the actual form of the derivative
908 {
909 number k_factorial = 2;
910 for (unsigned int k = 2; k <= n_derivatives; ++k)
911 {
912 for (unsigned int e = 0; e < n_entries; ++e)
913 values[k][e] *= k_factorial;
914 k_factorial *= static_cast<number>(k + 1);
915 }
916 }
917 break;
918
919 // manually implement case 0 (values only), case 1 (value + first
920 // derivative), and case 2 (up to second derivative) since they
921 // might be called often. then, we can unroll the inner loop and
922 // keep the temporary results as local variables to help the
923 // compiler with the pointer aliasing analysis.
924 case 0:
925 {
926 std::array<Number2, n_entries> value;
927 for (unsigned int e = 0; e < n_entries; ++e)
928 value[e] = weight;
929 for (unsigned int i = 0; i < n_supp; ++i)
930 for (unsigned int e = 0; e < n_entries; ++e)
931 value[e] *= (x[e] - lagrange_support_points[i]);
932
933 for (unsigned int e = 0; e < n_entries; ++e)
934 values[0][e] = value[e];
935 break;
936 }
937
938 case 1:
939 {
940 std::array<Number2, n_entries> value, derivative = {};
941 for (unsigned int e = 0; e < n_entries; ++e)
942 value[e] = weight;
943 for (unsigned int i = 0; i < n_supp; ++i)
944 for (unsigned int e = 0; e < n_entries; ++e)
945 {
946 const Number2 v = x[e] - lagrange_support_points[i];
947 derivative[e] = derivative[e] * v + value[e];
948 value[e] *= v;
949 }
950
951 for (unsigned int e = 0; e < n_entries; ++e)
952 {
953 values[0][e] = value[e];
954 values[1][e] = derivative[e];
955 }
956 break;
957 }
958
959 case 2:
960 {
961 std::array<Number2, n_entries> value, derivative = {},
962 second = {};
963 for (unsigned int e = 0; e < n_entries; ++e)
964 value[e] = weight;
965 for (unsigned int i = 0; i < n_supp; ++i)
966 for (unsigned int e = 0; e < n_entries; ++e)
967 {
968 const Number2 v = x[e] - lagrange_support_points[i];
969 second[e] = second[e] * v + derivative[e];
970 derivative[e] = derivative[e] * v + value[e];
971 value[e] *= v;
972 }
973
974 for (unsigned int e = 0; e < n_entries; ++e)
975 {
976 values[0][e] = value[e];
977 values[1][e] = derivative[e];
978 values[2][e] = static_cast<number>(2) * second[e];
979 }
980 break;
981 }
982 }
983 return;
984 }
985
986 Assert(coefficients.size() > 0, ExcEmptyObject());
987
988 // if derivatives are needed, then do it properly by the full
989 // Horner scheme
990 const unsigned int m = coefficients.size();
991 std::vector<std::array<Number2, n_entries>> a(coefficients.size());
992 for (unsigned int i = 0; i < coefficients.size(); ++i)
993 for (unsigned int e = 0; e < n_entries; ++e)
994 a[i][e] = coefficients[i];
995
996 unsigned int j_factorial = 1;
997
998 // loop over all requested derivatives. note that derivatives @p{j>m} are
999 // necessarily zero, as they differentiate the polynomial more often than
1000 // the highest power is
1001 const unsigned int min_valuessize_m = std::min(n_derivatives + 1, m);
1002 for (unsigned int j = 0; j < min_valuessize_m; ++j)
1003 {
1004 for (int k = m - 2; k >= static_cast<int>(j); --k)
1005 for (unsigned int e = 0; e < n_entries; ++e)
1006 a[k][e] += x[e] * a[k + 1][e];
1007 for (unsigned int e = 0; e < n_entries; ++e)
1008 values[j][e] = static_cast<number>(j_factorial) * a[j][e];
1009
1010 j_factorial *= j + 1;
1011 }
1012
1013 // fill higher derivatives by zero
1014 for (unsigned int j = min_valuessize_m; j <= n_derivatives; ++j)
1015 for (unsigned int e = 0; e < n_entries; ++e)
1016 values[j][e] = 0.;
1017 }
1018
1019
1020
1021 template <typename number>
1022 template <class Archive>
1023 inline void
1024 Polynomial<number>::serialize(Archive &ar, const unsigned int)
1025 {
1026 // forward to serialization function in the base class.
1027 ar &static_cast<Subscriptor &>(*this);
1028 ar &coefficients;
1029 ar &in_lagrange_product_form;
1030 ar &lagrange_support_points;
1031 ar &lagrange_weight;
1032 }
1033
1034
1035
1036 template <typename Number>
1037 Number
1038 jacobi_polynomial_value(const unsigned int degree,
1039 const int alpha,
1040 const int beta,
1041 const Number x)
1042 {
1043 Assert(alpha >= 0 && beta >= 0,
1044 ExcNotImplemented("Negative alpha/beta coefficients not supported"));
1045 // the Jacobi polynomial is evaluated using a recursion formula.
1046 Number p0, p1;
1047
1048 // The recursion formula is defined for the interval [-1, 1], so rescale
1049 // to that interval here
1050 const Number xeval = Number(-1) + 2. * x;
1051
1052 // initial values P_0(x), P_1(x):
1053 p0 = 1.0;
1054 if (degree == 0)
1055 return p0;
1056 p1 = ((alpha + beta + 2) * xeval + (alpha - beta)) / 2;
1057 if (degree == 1)
1058 return p1;
1059
1060 for (unsigned int i = 1; i < degree; ++i)
1061 {
1062 const Number v = 2 * i + (alpha + beta);
1063 const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
1064 const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
1065 const Number a3 = v * (v + 1) * (v + 2);
1066 const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
1067
1068 const Number pn = ((a2 + a3 * xeval) * p1 - a4 * p0) / a1;
1069 p0 = p1;
1070 p1 = pn;
1071 }
1072 return p1;
1073 }
1074
1075
1076
1077 template <typename Number>
1078 std::vector<Number>
1079 jacobi_polynomial_roots(const unsigned int degree,
1080 const int alpha,
1081 const int beta)
1082 {
1083 std::vector<Number> x(degree, 0.5);
1084
1085 // compute zeros with a Newton algorithm.
1086
1087 // Set tolerance. For long double we might not always get the additional
1088 // precision in a run time environment (e.g. with valgrind), so we must
1089 // limit the tolerance to double. Since we do a Newton iteration, doing
1090 // one more iteration after the residual has indicated convergence will be
1091 // enough for all number types due to the quadratic convergence of
1092 // Newton's method
1093
1094 const Number tolerance =
1095 4 * std::max(static_cast<Number>(std::numeric_limits<double>::epsilon()),
1096 std::numeric_limits<Number>::epsilon());
1097
1098 // The following implementation follows closely the one given in the
1099 // appendix of the book by Karniadakis and Sherwin: Spectral/hp element
1100 // methods for computational fluid dynamics (Oxford University Press,
1101 // 2005)
1102
1103 // If symmetric, we only need to compute the half of points
1104 const unsigned int n_points = (alpha == beta ? degree / 2 : degree);
1105 for (unsigned int k = 0; k < n_points; ++k)
1106 {
1107 // we take the zeros of the Chebyshev polynomial (alpha=beta=-0.5) as
1108 // initial values, corrected by the initial value
1109 Number r = 0.5 - 0.5 * std::cos(static_cast<Number>(2 * k + 1) /
1110 (2 * degree) * numbers::PI);
1111 if (k > 0)
1112 r = (r + x[k - 1]) / 2;
1113
1114 unsigned int converged = numbers::invalid_unsigned_int;
1115 for (unsigned int it = 1; it < 1000; ++it)
1116 {
1117 Number s = 0.;
1118 for (unsigned int i = 0; i < k; ++i)
1119 s += 1. / (r - x[i]);
1120
1121 // derivative of P_n^{alpha,beta}, rescaled to [0, 1]
1122 const Number J_x =
1123 (alpha + beta + degree + 1) *
1124 jacobi_polynomial_value(degree - 1, alpha + 1, beta + 1, r);
1125
1126 // value of P_n^{alpha,beta}
1127 const Number f = jacobi_polynomial_value(degree, alpha, beta, r);
1128 const Number delta = f / (f * s - J_x);
1129 r += delta;
1130 if (converged == numbers::invalid_unsigned_int &&
1131 std::abs(delta) < tolerance)
1132 converged = it;
1133
1134 // do one more iteration to ensure accuracy also for tighter
1135 // types than double (e.g. long double)
1136 if (it == converged + 1)
1137 break;
1138 }
1139
1141 ExcMessage("Newton iteration for zero of Jacobi polynomial "
1142 "did not converge."));
1143
1144 x[k] = r;
1145 }
1146
1147 // in case we assumed symmetry, fill up the missing values
1148 for (unsigned int k = n_points; k < degree; ++k)
1149 x[k] = 1.0 - x[degree - k - 1];
1150
1151 return x;
1152 }
1153
1154} // namespace Polynomials
1156
1157#endif
Definition point.h:111
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
static std::vector< std::unique_ptr< const std::vector< double > > > recursive_coefficients
Definition polynomial.h:574
static void compute_coefficients(const unsigned int p)
static const std::vector< double > & get_coefficients(const unsigned int p)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
static std::shared_mutex coefficients_lock
Definition polynomial.h:580
static void compute_coefficients(const unsigned int n, const unsigned int support_point, std::vector< double > &a)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
std::vector< double > compute_coefficients(const unsigned int p)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
static std::vector< Polynomial< number > > generate_complete_basis(const unsigned int degree)
static std::vector< number > make_vector(unsigned int n, const double coefficient)
number value(const number x) const
Definition polynomial.h:823
bool operator==(const Polynomial< number > &p) const
std::vector< number > coefficients
Definition polynomial.h:302
Polynomial< number > primitive() const
Polynomial< number > & operator+=(const Polynomial< number > &p)
void values_of_array(const std::array< Number2, n_entries > &points, const unsigned int n_derivatives, std::array< Number2, n_entries > *values) const
Definition polynomial.h:867
Polynomial< number > derivative() const
void transform_into_standard_form()
Definition polynomial.cc:96
void scale(const number factor)
Polynomial< number > & operator-=(const Polynomial< number > &p)
std::vector< number > lagrange_support_points
Definition polynomial.h:314
void shift(const number2 offset)
void print(std::ostream &out) const
void serialize(Archive &ar, const unsigned int version)
static void multiply(std::vector< number > &coefficients, const number factor)
void value(const Number2 x, const unsigned int n_derivatives, Number2 *values) const
Definition polynomial.h:853
Polynomial< number > & operator*=(const double s)
virtual std::size_t memory_consumption() const
unsigned int degree() const
Definition polynomial.h:806
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 2 > second
Definition grid_out.cc:4614
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcEmptyObject()
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
Number jacobi_polynomial_value(const unsigned int degree, const int alpha, const int beta, const Number x)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 > > &points)
std::vector< Number > jacobi_polynomial_roots(const unsigned int degree, const int alpha, const int beta)
static constexpr double PI
Definition numbers.h:259
static const unsigned int invalid_unsigned_int
Definition types.h:220
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)