Reference documentation for deal.II version GIT relicensing-761-gcc3619e707 2024-05-25 17:40: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#ifndef DEBUG
162#endif
163 void
164 values_of_array(const std::array<Number2, n_entries> &points,
165 const unsigned int n_derivatives,
166 std::array<Number2, n_entries> *values) const;
167
173 unsigned int
174 degree() const;
175
183 void
184 scale(const number factor);
185
201 template <typename number2>
202 void
203 shift(const number2 offset);
204
209 derivative() const;
210
216 primitive() const;
217
222 operator*=(const double s);
223
229
235
241
245 bool
246 operator==(const Polynomial<number> &p) const;
247
251 void
252 print(std::ostream &out) const;
253
259 template <class Archive>
260 void
261 serialize(Archive &ar, const unsigned int version);
262
266 virtual std::size_t
267 memory_consumption() const;
268
269 protected:
273 static void
274 scale(std::vector<number> &coefficients, const number factor);
275
279 template <typename number2>
280 static void
281 shift(std::vector<number> &coefficients, const number2 shift);
282
286 static void
287 multiply(std::vector<number> &coefficients, const number factor);
288
294 void
296
305 std::vector<number> coefficients;
306
312
317 std::vector<number> lagrange_support_points;
318
324 };
325
326
331 template <typename number>
332 class Monomial : public Polynomial<number>
333 {
334 public:
339 Monomial(const unsigned int n, const double coefficient = 1.);
340
347 static std::vector<Polynomial<number>>
348 generate_complete_basis(const unsigned int degree);
349
350 private:
354 static std::vector<number>
355 make_vector(unsigned int n, const double coefficient);
356 };
357
358
374 class LagrangeEquidistant : public Polynomial<double>
375 {
376 public:
382 LagrangeEquidistant(const unsigned int n, const unsigned int support_point);
383
392 static std::vector<Polynomial<double>>
393 generate_complete_basis(const unsigned int degree);
394
395 private:
400 static void
401 compute_coefficients(const unsigned int n,
402 const unsigned int support_point,
403 std::vector<double> &a);
404 };
405
406
407
414 std::vector<Polynomial<double>>
415 generate_complete_Lagrange_basis(const std::vector<Point<1>> &points);
416
417
418
432 class Legendre : public Polynomial<double>
433 {
434 public:
438 Legendre(const unsigned int p);
439
446 static std::vector<Polynomial<double>>
447 generate_complete_basis(const unsigned int degree);
448 };
449
469 class Lobatto : public Polynomial<double>
470 {
471 public:
476 Lobatto(const unsigned int p = 0);
477
482 static std::vector<Polynomial<double>>
483 generate_complete_basis(const unsigned int p);
484
485 private:
489 std::vector<double>
490 compute_coefficients(const unsigned int p);
491 };
492
493
494
532 class Hierarchical : public Polynomial<double>
533 {
534 public:
539 Hierarchical(const unsigned int p);
540
551 static std::vector<Polynomial<double>>
552 generate_complete_basis(const unsigned int degree);
553
554 private:
558 static void
559 compute_coefficients(const unsigned int p);
560
565 static const std::vector<double> &
566 get_coefficients(const unsigned int p);
567
576 static std::vector<std::unique_ptr<const std::vector<double>>>
578
583 static std::shared_mutex coefficients_lock;
584 };
585
586
587
615 class HermiteInterpolation : public Polynomial<double>
616 {
617 public:
622 HermiteInterpolation(const unsigned int p);
623
629 static std::vector<Polynomial<double>>
630 generate_complete_basis(const unsigned int p);
631 };
632
633
634
737 {
738 public:
743 HermiteLikeInterpolation(const unsigned int degree,
744 const unsigned int index);
745
750 static std::vector<Polynomial<double>>
751 generate_complete_basis(const unsigned int degree);
752 };
753
754
755
756 /*
757 * Evaluate a Jacobi polynomial @f$ P_n^{\alpha, \beta}(x) @f$ specified by the
758 * parameters @p alpha, @p beta, @p n, where @p n is the degree of the
759 * Jacobi polynomial.
760 *
761 * @note The Jacobi polynomials are not orthonormal and are defined on the
762 * unit interval @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often
763 * used in literature. @p x is the point of evaluation.
764 */
765 template <typename Number>
766 Number
767 jacobi_polynomial_value(const unsigned int degree,
768 const int alpha,
769 const int beta,
770 const Number x);
771
772
785 template <typename Number>
786 std::vector<Number>
787 jacobi_polynomial_roots(const unsigned int degree,
788 const int alpha,
789 const int beta);
790} // namespace Polynomials
791
792
795/* -------------------------- inline functions --------------------- */
796
797namespace Polynomials
798{
799 template <typename number>
801 : in_lagrange_product_form(false)
802 , lagrange_weight(1.)
803 {}
804
805
806
807 template <typename number>
808 inline unsigned int
810 {
811 if (in_lagrange_product_form == true)
812 {
813 return lagrange_support_points.size();
814 }
815 else
816 {
817 Assert(coefficients.size() > 0, ExcEmptyObject());
818 return coefficients.size() - 1;
819 }
820 }
821
822
823
824 template <typename number>
825 inline number
826 Polynomial<number>::value(const number x) const
827 {
828 if (in_lagrange_product_form == false)
829 {
830 Assert(coefficients.size() > 0, ExcEmptyObject());
831
832 // Horner scheme
833 const unsigned int m = coefficients.size();
834 number value = coefficients.back();
835 for (int k = m - 2; k >= 0; --k)
836 value = value * x + coefficients[k];
837 return value;
838 }
839 else
840 {
841 // direct evaluation of Lagrange polynomial
842 const unsigned int m = lagrange_support_points.size();
843 number value = 1.;
844 for (unsigned int j = 0; j < m; ++j)
845 value *= x - lagrange_support_points[j];
846 value *= lagrange_weight;
847 return value;
848 }
849 }
850
851
852
853 template <typename number>
854 template <typename Number2>
855 inline void
857 const unsigned int n_derivatives,
858 Number2 *values) const
859 {
860 values_of_array(std::array<Number2, 1ul>{{x}},
861 n_derivatives,
862 reinterpret_cast<std::array<Number2, 1ul> *>(values));
863 }
864
865
866
867 template <typename number>
868 template <std::size_t n_entries, typename Number2>
869 inline
870#ifndef DEBUG
872#endif
873 void
875 const std::array<Number2, n_entries> &x,
876 const unsigned int n_derivatives,
877 std::array<Number2, n_entries> *values) const
878 {
879 // evaluate Lagrange polynomial and derivatives
880 if (in_lagrange_product_form == true)
881 {
882 // to compute the value and all derivatives of a polynomial of the
883 // form (x-x_1)*(x-x_2)*...*(x-x_n), expand the derivatives like
884 // automatic differentiation does.
885 const unsigned int n_supp = lagrange_support_points.size();
886 const number weight = lagrange_weight;
887 switch (n_derivatives)
888 {
889 default:
890 for (unsigned int e = 0; e < n_entries; ++e)
891 values[0][e] = weight;
892 for (unsigned int k = 1; k <= n_derivatives; ++k)
893 for (unsigned int e = 0; e < n_entries; ++e)
894 values[k][e] = 0.;
895 for (unsigned int i = 0; i < n_supp; ++i)
896 {
897 std::array<Number2, n_entries> v = x;
898 for (unsigned int e = 0; e < n_entries; ++e)
899 v[e] -= lagrange_support_points[i];
900
901 // multiply by (x-x_i) and compute action on all derivatives,
902 // too (inspired from automatic differentiation: implement the
903 // product rule for the old value and the new variable 'v',
904 // i.e., expand value v and derivative one). since we reuse a
905 // value from the next lower derivative from the steps before,
906 // need to start from the highest derivative
907 for (unsigned int k = n_derivatives; k > 0; --k)
908 for (unsigned int e = 0; e < n_entries; ++e)
909 values[k][e] = (values[k][e] * v[e] + values[k - 1][e]);
910 for (unsigned int e = 0; e < n_entries; ++e)
911 values[0][e] *= v[e];
912 }
913 // finally, multiply derivatives by k! to transform the product
914 // p_n = p^(n)(x)/k! into the actual form of the derivative
915 {
916 number k_factorial = 2;
917 for (unsigned int k = 2; k <= n_derivatives; ++k)
918 {
919 for (unsigned int e = 0; e < n_entries; ++e)
920 values[k][e] *= k_factorial;
921 k_factorial *= static_cast<number>(k + 1);
922 }
923 }
924 break;
925
926 // manually implement case 0 (values only), case 1 (value + first
927 // derivative), and case 2 (up to second derivative) since they
928 // might be called often. then, we can unroll the inner loop and
929 // keep the temporary results as local variables to help the
930 // compiler with the pointer aliasing analysis.
931 case 0:
932 {
933 std::array<Number2, n_entries> value;
934 for (unsigned int e = 0; e < n_entries; ++e)
935 value[e] = weight;
936 for (unsigned int i = 0; i < n_supp; ++i)
937 for (unsigned int e = 0; e < n_entries; ++e)
938 value[e] *= (x[e] - lagrange_support_points[i]);
939
940 for (unsigned int e = 0; e < n_entries; ++e)
941 values[0][e] = value[e];
942 break;
943 }
944
945 case 1:
946 {
947 std::array<Number2, n_entries> value, derivative = {};
948 for (unsigned int e = 0; e < n_entries; ++e)
949 value[e] = weight;
950 for (unsigned int i = 0; i < n_supp; ++i)
951 for (unsigned int e = 0; e < n_entries; ++e)
952 {
953 const Number2 v = x[e] - lagrange_support_points[i];
954 derivative[e] = derivative[e] * v + value[e];
955 value[e] *= v;
956 }
957
958 for (unsigned int e = 0; e < n_entries; ++e)
959 {
960 values[0][e] = value[e];
961 values[1][e] = derivative[e];
962 }
963 break;
964 }
965
966 case 2:
967 {
968 std::array<Number2, n_entries> value, derivative = {},
969 second = {};
970 for (unsigned int e = 0; e < n_entries; ++e)
971 value[e] = weight;
972 for (unsigned int i = 0; i < n_supp; ++i)
973 for (unsigned int e = 0; e < n_entries; ++e)
974 {
975 const Number2 v = x[e] - lagrange_support_points[i];
976 second[e] = second[e] * v + derivative[e];
977 derivative[e] = derivative[e] * v + value[e];
978 value[e] *= v;
979 }
980
981 for (unsigned int e = 0; e < n_entries; ++e)
982 {
983 values[0][e] = value[e];
984 values[1][e] = derivative[e];
985 values[2][e] = static_cast<number>(2) * second[e];
986 }
987 break;
988 }
989 }
990 return;
991 }
992
993 Assert(coefficients.size() > 0, ExcEmptyObject());
994
995 // if derivatives are needed, then do it properly by the full
996 // Horner scheme
997 const unsigned int m = coefficients.size();
998 std::vector<std::array<Number2, n_entries>> a(coefficients.size());
999 for (unsigned int i = 0; i < coefficients.size(); ++i)
1000 for (unsigned int e = 0; e < n_entries; ++e)
1001 a[i][e] = coefficients[i];
1002
1003 unsigned int j_factorial = 1;
1004
1005 // loop over all requested derivatives. note that derivatives @p{j>m} are
1006 // necessarily zero, as they differentiate the polynomial more often than
1007 // the highest power is
1008 const unsigned int min_valuessize_m = std::min(n_derivatives + 1, m);
1009 for (unsigned int j = 0; j < min_valuessize_m; ++j)
1010 {
1011 for (int k = m - 2; k >= static_cast<int>(j); --k)
1012 for (unsigned int e = 0; e < n_entries; ++e)
1013 a[k][e] += x[e] * a[k + 1][e];
1014 for (unsigned int e = 0; e < n_entries; ++e)
1015 values[j][e] = static_cast<number>(j_factorial) * a[j][e];
1016
1017 j_factorial *= j + 1;
1018 }
1019
1020 // fill higher derivatives by zero
1021 for (unsigned int j = min_valuessize_m; j <= n_derivatives; ++j)
1022 for (unsigned int e = 0; e < n_entries; ++e)
1023 values[j][e] = 0.;
1024 }
1025
1026
1027
1028 template <typename number>
1029 template <class Archive>
1030 inline void
1031 Polynomial<number>::serialize(Archive &ar, const unsigned int)
1032 {
1033 // forward to serialization function in the base class.
1034 ar &static_cast<Subscriptor &>(*this);
1035 ar &coefficients;
1036 ar &in_lagrange_product_form;
1037 ar &lagrange_support_points;
1038 ar &lagrange_weight;
1039 }
1040
1041
1042
1043 template <typename Number>
1044 Number
1045 jacobi_polynomial_value(const unsigned int degree,
1046 const int alpha,
1047 const int beta,
1048 const Number x)
1049 {
1050 Assert(alpha >= 0 && beta >= 0,
1051 ExcNotImplemented("Negative alpha/beta coefficients not supported"));
1052 // the Jacobi polynomial is evaluated using a recursion formula.
1053 Number p0, p1;
1054
1055 // The recursion formula is defined for the interval [-1, 1], so rescale
1056 // to that interval here
1057 const Number xeval = Number(-1) + 2. * x;
1058
1059 // initial values P_0(x), P_1(x):
1060 p0 = 1.0;
1061 if (degree == 0)
1062 return p0;
1063 p1 = ((alpha + beta + 2) * xeval + (alpha - beta)) / 2;
1064 if (degree == 1)
1065 return p1;
1066
1067 for (unsigned int i = 1; i < degree; ++i)
1068 {
1069 const Number v = 2 * i + (alpha + beta);
1070 const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
1071 const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
1072 const Number a3 = v * (v + 1) * (v + 2);
1073 const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
1074
1075 const Number pn = ((a2 + a3 * xeval) * p1 - a4 * p0) / a1;
1076 p0 = p1;
1077 p1 = pn;
1078 }
1079 return p1;
1080 }
1081
1082
1083
1084 template <typename Number>
1085 std::vector<Number>
1086 jacobi_polynomial_roots(const unsigned int degree,
1087 const int alpha,
1088 const int beta)
1089 {
1090 std::vector<Number> x(degree, 0.5);
1091
1092 // compute zeros with a Newton algorithm.
1093
1094 // Set tolerance. For long double we might not always get the additional
1095 // precision in a run time environment (e.g. with valgrind), so we must
1096 // limit the tolerance to double. Since we do a Newton iteration, doing
1097 // one more iteration after the residual has indicated convergence will be
1098 // enough for all number types due to the quadratic convergence of
1099 // Newton's method
1100
1101 const Number tolerance =
1102 4 * std::max(static_cast<Number>(std::numeric_limits<double>::epsilon()),
1103 std::numeric_limits<Number>::epsilon());
1104
1105 // The following implementation follows closely the one given in the
1106 // appendix of the book by Karniadakis and Sherwin: Spectral/hp element
1107 // methods for computational fluid dynamics (Oxford University Press,
1108 // 2005)
1109
1110 // If symmetric, we only need to compute the half of points
1111 const unsigned int n_points = (alpha == beta ? degree / 2 : degree);
1112 for (unsigned int k = 0; k < n_points; ++k)
1113 {
1114 // we take the zeros of the Chebyshev polynomial (alpha=beta=-0.5) as
1115 // initial values, corrected by the initial value
1116 Number r = 0.5 - 0.5 * std::cos(static_cast<Number>(2 * k + 1) /
1117 (2 * degree) * numbers::PI);
1118 if (k > 0)
1119 r = (r + x[k - 1]) / 2;
1120
1121 unsigned int converged = numbers::invalid_unsigned_int;
1122 for (unsigned int it = 1; it < 1000; ++it)
1123 {
1124 Number s = 0.;
1125 for (unsigned int i = 0; i < k; ++i)
1126 s += 1. / (r - x[i]);
1127
1128 // derivative of P_n^{alpha,beta}, rescaled to [0, 1]
1129 const Number J_x =
1130 (alpha + beta + degree + 1) *
1131 jacobi_polynomial_value(degree - 1, alpha + 1, beta + 1, r);
1132
1133 // value of P_n^{alpha,beta}
1134 const Number f = jacobi_polynomial_value(degree, alpha, beta, r);
1135 const Number delta = f / (f * s - J_x);
1136 r += delta;
1137 if (converged == numbers::invalid_unsigned_int &&
1138 std::abs(delta) < tolerance)
1139 converged = it;
1140
1141 // do one more iteration to ensure accuracy also for tighter
1142 // types than double (e.g. long double)
1143 if (it == converged + 1)
1144 break;
1145 }
1146
1148 ExcMessage("Newton iteration for zero of Jacobi polynomial "
1149 "did not converge."));
1150
1151 x[k] = r;
1152 }
1153
1154 // in case we assumed symmetry, fill up the missing values
1155 for (unsigned int k = n_points; k < degree; ++k)
1156 x[k] = 1.0 - x[degree - k - 1];
1157
1158 return x;
1159 }
1160
1161} // namespace Polynomials
1163
1164#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:577
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:583
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:826
bool operator==(const Polynomial< number > &p) const
std::vector< number > coefficients
Definition polynomial.h:305
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:874
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:317
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:856
Polynomial< number > & operator*=(const double s)
virtual std::size_t memory_consumption() const
unsigned int degree() const
Definition polynomial.h:809
#define DEAL_II_ALWAYS_INLINE
Definition config.h:109
#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:4624
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 > &)