Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3075-gc235bd4825 2025-04-15 08:10:00+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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
polynomial.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2000 - 2024 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
25#include <deal.II/base/point.h>
26#include <deal.II/base/types.h>
27
28#include <array>
29#include <limits>
30#include <memory>
31#include <shared_mutex>
32#include <vector>
33
35
45namespace Polynomials
46{
66 template <typename number>
68 {
69 public:
78 Polynomial(const std::vector<number> &coefficients);
79
83 Polynomial(const unsigned int n);
84
93 const unsigned int evaluation_point);
94
99
109 number
110 value(const number x) const;
111
122 void
123 value(const number x, std::vector<number> &values) const;
124
143 template <typename Number2>
144 void
145 value(const Number2 x,
146 const unsigned int n_derivatives,
147 Number2 *values) const;
148
161 template <std::size_t n_entries, typename Number2>
162#ifndef DEBUG
164#endif
165 void
166 values_of_array(const std::array<Number2, n_entries> &points,
167 const unsigned int n_derivatives,
168 std::array<Number2, n_entries> *values) const;
169
175 unsigned int
176 degree() const;
177
185 void
186 scale(const number factor);
187
203 template <typename number2>
204 void
205 shift(const number2 offset);
206
211 derivative() const;
212
218 primitive() const;
219
224 operator*=(const double s);
225
231
237
243
247 bool
248 operator==(const Polynomial<number> &p) const;
249
253 void
254 print(std::ostream &out) const;
255
261 template <class Archive>
262 void
263 serialize(Archive &ar, const unsigned int version);
264
268 virtual std::size_t
269 memory_consumption() const;
270
271 protected:
275 static void
276 scale(std::vector<number> &coefficients, const number factor);
277
281 template <typename number2>
282 static void
283 shift(std::vector<number> &coefficients, const number2 shift);
284
288 static void
289 multiply(std::vector<number> &coefficients, const number factor);
290
296 void
298
307 std::vector<number> coefficients;
308
314
319 std::vector<number> lagrange_support_points;
320
326 };
327
328
333 template <typename number>
334 class Monomial : public Polynomial<number>
335 {
336 public:
341 Monomial(const unsigned int n, const double coefficient = 1.);
342
349 static std::vector<Polynomial<number>>
350 generate_complete_basis(const unsigned int degree);
351
352 private:
356 static std::vector<number>
357 make_vector(unsigned int n, const double coefficient);
358 };
359
360
376 class LagrangeEquidistant : public Polynomial<double>
377 {
378 public:
384 LagrangeEquidistant(const unsigned int n, const unsigned int support_point);
385
394 static std::vector<Polynomial<double>>
395 generate_complete_basis(const unsigned int degree);
396
397 private:
402 static void
403 compute_coefficients(const unsigned int n,
404 const unsigned int support_point,
405 std::vector<double> &a);
406 };
407
408
409
416 std::vector<Polynomial<double>>
417 generate_complete_Lagrange_basis(const std::vector<Point<1>> &points);
418
419
420
434 class Legendre : public Polynomial<double>
435 {
436 public:
440 Legendre(const unsigned int p);
441
448 static std::vector<Polynomial<double>>
449 generate_complete_basis(const unsigned int degree);
450 };
451
471 class Lobatto : public Polynomial<double>
472 {
473 public:
478 Lobatto(const unsigned int p = 0);
479
484 static std::vector<Polynomial<double>>
485 generate_complete_basis(const unsigned int p);
486
487 private:
491 std::vector<double>
492 compute_coefficients(const unsigned int p);
493 };
494
495
496
534 class Hierarchical : public Polynomial<double>
535 {
536 public:
541 Hierarchical(const unsigned int p);
542
553 static std::vector<Polynomial<double>>
554 generate_complete_basis(const unsigned int degree);
555
556 private:
560 static void
561 compute_coefficients(const unsigned int p);
562
567 static const std::vector<double> &
568 get_coefficients(const unsigned int p);
569
578 static std::vector<std::unique_ptr<const std::vector<double>>>
580
585 static std::shared_mutex coefficients_lock;
586 };
587
588
589
617 class HermiteInterpolation : public Polynomial<double>
618 {
619 public:
624 HermiteInterpolation(const unsigned int p);
625
631 static std::vector<Polynomial<double>>
632 generate_complete_basis(const unsigned int p);
633 };
634
635
636
739 {
740 public:
745 HermiteLikeInterpolation(const unsigned int degree,
746 const unsigned int index);
747
752 static std::vector<Polynomial<double>>
753 generate_complete_basis(const unsigned int degree);
754 };
755
756
757
758 /*
759 * Evaluate a Jacobi polynomial @f$ P_n^{\alpha, \beta}(x) @f$ specified by the
760 * parameters @p alpha, @p beta, @p n, where @p n is the degree of the
761 * Jacobi polynomial.
762 *
763 * @note The Jacobi polynomials are not orthonormal and are defined on the
764 * unit interval @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often
765 * used in literature. @p x is the point of evaluation.
766 */
767 template <typename Number>
768 Number
769 jacobi_polynomial_value(const unsigned int degree,
770 const int alpha,
771 const int beta,
772 const Number x);
773
774
787 template <typename Number>
788 std::vector<Number>
789 jacobi_polynomial_roots(const unsigned int degree,
790 const int alpha,
791 const int beta);
792} // namespace Polynomials
793
794
797/* -------------------------- inline functions --------------------- */
798
799namespace Polynomials
800{
801 template <typename number>
803 : in_lagrange_product_form(false)
804 , lagrange_weight(1.)
805 {}
806
807
808
809 template <typename number>
810 inline unsigned int
812 {
813 if (in_lagrange_product_form == true)
814 {
815 return lagrange_support_points.size();
816 }
817 else
818 {
819 Assert(coefficients.size() > 0, ExcEmptyObject());
820 return coefficients.size() - 1;
821 }
822 }
823
824
825
826 template <typename number>
827 inline number
828 Polynomial<number>::value(const number x) const
829 {
830 if (in_lagrange_product_form == false)
831 {
832 Assert(coefficients.size() > 0, ExcEmptyObject());
833
834 // Horner scheme
835 const unsigned int m = coefficients.size();
836 number value = coefficients.back();
837 for (int k = m - 2; k >= 0; --k)
838 value = value * x + coefficients[k];
839 return value;
840 }
841 else
842 {
843 // direct evaluation of Lagrange polynomial
844 const unsigned int m = lagrange_support_points.size();
845 number value = 1.;
846 for (unsigned int j = 0; j < m; ++j)
847 value *= x - lagrange_support_points[j];
848 value *= lagrange_weight;
849 return value;
850 }
851 }
852
853
854
855 template <typename number>
856 template <typename Number2>
857 inline void
859 const unsigned int n_derivatives,
860 Number2 *values) const
861 {
862 values_of_array(std::array<Number2, 1ul>{{x}},
863 n_derivatives,
864 reinterpret_cast<std::array<Number2, 1ul> *>(values));
865 }
866
867
868
869 template <typename number>
870 template <std::size_t n_entries, typename Number2>
871 inline
872#ifndef DEBUG
874#endif
875 void
877 const std::array<Number2, n_entries> &x,
878 const unsigned int n_derivatives,
879 std::array<Number2, n_entries> *values) const
880 {
881 // evaluate Lagrange polynomial and derivatives
882 if (in_lagrange_product_form == true)
883 {
884 // to compute the value and all derivatives of a polynomial of the
885 // form (x-x_1)*(x-x_2)*...*(x-x_n), expand the derivatives like
886 // automatic differentiation does.
887 const unsigned int n_supp = lagrange_support_points.size();
888 const number weight = lagrange_weight;
889 switch (n_derivatives)
890 {
891 default:
892 for (unsigned int e = 0; e < n_entries; ++e)
893 values[0][e] = weight;
894 for (unsigned int k = 1; k <= n_derivatives; ++k)
895 for (unsigned int e = 0; e < n_entries; ++e)
896 values[k][e] = 0.;
897 for (unsigned int i = 0; i < n_supp; ++i)
898 {
899 std::array<Number2, n_entries> v = x;
900 for (unsigned int e = 0; e < n_entries; ++e)
901 v[e] -= lagrange_support_points[i];
902
903 // multiply by (x-x_i) and compute action on all derivatives,
904 // too (inspired from automatic differentiation: implement the
905 // product rule for the old value and the new variable 'v',
906 // i.e., expand value v and derivative one). since we reuse a
907 // value from the next lower derivative from the steps before,
908 // need to start from the highest derivative
909 for (unsigned int k = n_derivatives; k > 0; --k)
910 for (unsigned int e = 0; e < n_entries; ++e)
911 values[k][e] = (values[k][e] * v[e] + values[k - 1][e]);
912 for (unsigned int e = 0; e < n_entries; ++e)
913 values[0][e] *= v[e];
914 }
915 // finally, multiply derivatives by k! to transform the product
916 // p_n = p^(n)(x)/k! into the actual form of the derivative
917 {
918 number k_factorial = 2;
919 for (unsigned int k = 2; k <= n_derivatives; ++k)
920 {
921 for (unsigned int e = 0; e < n_entries; ++e)
922 values[k][e] *= k_factorial;
923 k_factorial *= static_cast<number>(k + 1);
924 }
925 }
926 break;
927
928 // manually implement case 0 (values only), case 1 (value + first
929 // derivative), and case 2 (up to second derivative) since they
930 // might be called often. then, we can unroll the inner loop and
931 // keep the temporary results as local variables to help the
932 // compiler with the pointer aliasing analysis.
933 case 0:
934 {
935 std::array<Number2, n_entries> value;
936 for (unsigned int e = 0; e < n_entries; ++e)
937 value[e] = weight;
938 for (unsigned int i = 0; i < n_supp; ++i)
939 for (unsigned int e = 0; e < n_entries; ++e)
940 value[e] *= (x[e] - lagrange_support_points[i]);
941
942 for (unsigned int e = 0; e < n_entries; ++e)
943 values[0][e] = value[e];
944 break;
945 }
946
947 case 1:
948 {
949 std::array<Number2, n_entries> value, derivative = {};
950 for (unsigned int e = 0; e < n_entries; ++e)
951 value[e] = weight;
952 for (unsigned int i = 0; i < n_supp; ++i)
953 for (unsigned int e = 0; e < n_entries; ++e)
954 {
955 const Number2 v = x[e] - lagrange_support_points[i];
956 derivative[e] = derivative[e] * v + value[e];
957 value[e] *= v;
958 }
959
960 for (unsigned int e = 0; e < n_entries; ++e)
961 {
962 values[0][e] = value[e];
963 values[1][e] = derivative[e];
964 }
965 break;
966 }
967
968 case 2:
969 {
970 std::array<Number2, n_entries> value, derivative = {},
971 second = {};
972 for (unsigned int e = 0; e < n_entries; ++e)
973 value[e] = weight;
974 for (unsigned int i = 0; i < n_supp; ++i)
975 for (unsigned int e = 0; e < n_entries; ++e)
976 {
977 const Number2 v = x[e] - lagrange_support_points[i];
978 second[e] = second[e] * v + derivative[e];
979 derivative[e] = derivative[e] * v + value[e];
980 value[e] *= v;
981 }
982
983 for (unsigned int e = 0; e < n_entries; ++e)
984 {
985 values[0][e] = value[e];
986 values[1][e] = derivative[e];
987 values[2][e] = static_cast<number>(2) * second[e];
988 }
989 break;
990 }
991 }
992 return;
993 }
994
995 Assert(coefficients.size() > 0, ExcEmptyObject());
996
997 // if derivatives are needed, then do it properly by the full
998 // Horner scheme
999 const unsigned int m = coefficients.size();
1000 std::vector<std::array<Number2, n_entries>> a(coefficients.size());
1001 for (unsigned int i = 0; i < coefficients.size(); ++i)
1002 for (unsigned int e = 0; e < n_entries; ++e)
1003 a[i][e] = coefficients[i];
1004
1005 unsigned int j_factorial = 1;
1006
1007 // loop over all requested derivatives. note that derivatives @p{j>m} are
1008 // necessarily zero, as they differentiate the polynomial more often than
1009 // the highest power is
1010 const unsigned int min_valuessize_m = std::min(n_derivatives + 1, m);
1011 for (unsigned int j = 0; j < min_valuessize_m; ++j)
1012 {
1013 for (int k = m - 2; k >= static_cast<int>(j); --k)
1014 for (unsigned int e = 0; e < n_entries; ++e)
1015 a[k][e] += x[e] * a[k + 1][e];
1016 for (unsigned int e = 0; e < n_entries; ++e)
1017 values[j][e] = static_cast<number>(j_factorial) * a[j][e];
1018
1019 j_factorial *= j + 1;
1020 }
1021
1022 // fill higher derivatives by zero
1023 for (unsigned int j = min_valuessize_m; j <= n_derivatives; ++j)
1024 for (unsigned int e = 0; e < n_entries; ++e)
1025 values[j][e] = 0.;
1026 }
1027
1028
1029
1030 template <typename number>
1031 template <class Archive>
1032 inline void
1033 Polynomial<number>::serialize(Archive &ar, const unsigned int)
1034 {
1035 // forward to serialization function in the base class.
1036 ar &static_cast<EnableObserverPointer &>(*this);
1037 ar &coefficients;
1038 ar &in_lagrange_product_form;
1039 ar &lagrange_support_points;
1040 ar &lagrange_weight;
1041 }
1042
1043
1044
1045 template <typename Number>
1046 Number
1047 jacobi_polynomial_value(const unsigned int degree,
1048 const int alpha,
1049 const int beta,
1050 const Number x)
1051 {
1052 Assert(alpha >= 0 && beta >= 0,
1053 ExcNotImplemented("Negative alpha/beta coefficients not supported"));
1054 // the Jacobi polynomial is evaluated using a recursion formula.
1055 Number p0, p1;
1056
1057 // The recursion formula is defined for the interval [-1, 1], so rescale
1058 // to that interval here
1059 const Number xeval = Number(-1) + 2. * x;
1060
1061 // initial values P_0(x), P_1(x):
1062 p0 = 1.0;
1063 if (degree == 0)
1064 return p0;
1065 p1 = ((alpha + beta + 2) * xeval + (alpha - beta)) / 2;
1066 if (degree == 1)
1067 return p1;
1068
1069 for (unsigned int i = 1; i < degree; ++i)
1070 {
1071 const Number v = 2 * i + (alpha + beta);
1072 const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
1073 const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
1074 const Number a3 = v * (v + 1) * (v + 2);
1075 const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
1076
1077 const Number pn = ((a2 + a3 * xeval) * p1 - a4 * p0) / a1;
1078 p0 = p1;
1079 p1 = pn;
1080 }
1081 return p1;
1082 }
1083
1084
1085
1086 template <typename Number>
1087 std::vector<Number>
1088 jacobi_polynomial_roots(const unsigned int degree,
1089 const int alpha,
1090 const int beta)
1091 {
1092 std::vector<Number> x(degree, 0.5);
1093
1094 // compute zeros with a Newton algorithm.
1095
1096 // Set tolerance. For long double we might not always get the additional
1097 // precision in a run time environment (e.g. with valgrind), so we must
1098 // limit the tolerance to double. Since we do a Newton iteration, doing
1099 // one more iteration after the residual has indicated convergence will be
1100 // enough for all number types due to the quadratic convergence of
1101 // Newton's method
1102
1103 const Number tolerance =
1104 4 * std::max(static_cast<Number>(std::numeric_limits<double>::epsilon()),
1105 std::numeric_limits<Number>::epsilon());
1106
1107 // The following implementation follows closely the one given in the
1108 // appendix of the book by Karniadakis and Sherwin: Spectral/hp element
1109 // methods for computational fluid dynamics (Oxford University Press,
1110 // 2005)
1111
1112 // If symmetric, we only need to compute the half of points
1113 const unsigned int n_points = (alpha == beta ? degree / 2 : degree);
1114 for (unsigned int k = 0; k < n_points; ++k)
1115 {
1116 // we take the zeros of the Chebyshev polynomial (alpha=beta=-0.5) as
1117 // initial values, corrected by the initial value
1118 Number r = 0.5 - 0.5 * std::cos(static_cast<Number>(2 * k + 1) /
1119 (2 * degree) * numbers::PI);
1120 if (k > 0)
1121 r = (r + x[k - 1]) / 2;
1122
1123 unsigned int converged = numbers::invalid_unsigned_int;
1124 for (unsigned int it = 1; it < 1000; ++it)
1125 {
1126 Number s = 0.;
1127 for (unsigned int i = 0; i < k; ++i)
1128 s += 1. / (r - x[i]);
1129
1130 // derivative of P_n^{alpha,beta}, rescaled to [0, 1]
1131 const Number J_x =
1132 (alpha + beta + degree + 1) *
1133 jacobi_polynomial_value(degree - 1, alpha + 1, beta + 1, r);
1134
1135 // value of P_n^{alpha,beta}
1136 const Number f = jacobi_polynomial_value(degree, alpha, beta, r);
1137 const Number delta = f / (f * s - J_x);
1138 r += delta;
1139 if (converged == numbers::invalid_unsigned_int &&
1140 std::abs(delta) < tolerance)
1141 converged = it;
1142
1143 // do one more iteration to ensure accuracy also for tighter
1144 // types than double (e.g. long double)
1145 if (it == converged + 1)
1146 break;
1147 }
1148
1150 ExcMessage("Newton iteration for zero of Jacobi polynomial "
1151 "did not converge."));
1152
1153 x[k] = r;
1154 }
1155
1156 // in case we assumed symmetry, fill up the missing values
1157 for (unsigned int k = n_points; k < degree; ++k)
1158 x[k] = 1.0 - x[degree - k - 1];
1159
1160 return x;
1161 }
1162
1163} // namespace Polynomials
1165
1166#endif
Definition point.h:113
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:579
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:585
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:828
bool operator==(const Polynomial< number > &p) const
std::vector< number > coefficients
Definition polynomial.h:307
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:876
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:319
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:858
Polynomial< number > & operator*=(const double s)
virtual std::size_t memory_consumption() const
unsigned int degree() const
Definition polynomial.h:811
#define DEAL_II_ALWAYS_INLINE
Definition config.h:161
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
Point< 2 > second
Definition grid_out.cc:4633
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)
constexpr double PI
Definition numbers.h:240
constexpr unsigned int invalid_unsigned_int
Definition types.h:238
::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 > &)