16 #ifndef dealii_polynomial_h 17 #define dealii_polynomial_h 59 template <
typename number>
86 const unsigned int evaluation_point);
103 value(
const number x)
const;
116 value(
const number x, std::vector<number> &
values)
const;
136 template <
typename Number2>
138 value(
const Number2 x,
139 const unsigned int n_derivatives,
140 Number2 * values)
const;
158 scale(
const number factor);
175 template <
typename number2>
177 shift(
const number2 offset);
226 print(std::ostream &out)
const;
233 template <
class Archive>
235 serialize(Archive &ar,
const unsigned int version);
248 scale(std::vector<number> &coefficients,
const number factor);
253 template <
typename number2>
255 shift(std::vector<number> &coefficients,
const number2 shift);
261 multiply(std::vector<number> &coefficients,
const number factor);
305 template <
typename number>
313 Monomial(
const unsigned int n,
const double coefficient = 1.);
321 static std::vector<Polynomial<number>>
322 generate_complete_basis(
const unsigned int degree);
328 static std::vector<number>
329 make_vector(
unsigned int n,
const double coefficient);
366 static std::vector<Polynomial<double>>
367 generate_complete_basis(
const unsigned int degree);
375 compute_coefficients(
const unsigned int n,
376 const unsigned int support_point,
377 std::vector<double> &a);
388 std::vector<Polynomial<double>>
420 static std::vector<Polynomial<double>>
421 generate_complete_basis(
const unsigned int degree);
450 Lobatto(
const unsigned int p = 0);
456 static std::vector<Polynomial<double>>
457 generate_complete_basis(
const unsigned int p);
464 compute_coefficients(
const unsigned int p);
525 static std::vector<Polynomial<double>>
526 generate_complete_basis(
const unsigned int degree);
533 compute_coefficients(
const unsigned int p);
539 static const std::vector<double> &
540 get_coefficients(
const unsigned int p);
550 static std::vector<std::unique_ptr<const std::vector<double>>>
597 static std::vector<Polynomial<double>>
598 generate_complete_basis(
const unsigned int p);
712 const unsigned int index);
718 static std::vector<Polynomial<double>>
719 generate_complete_basis(
const unsigned int degree);
733 template <
typename Number>
753 template <
typename Number>
767 template <
typename number>
775 template <
typename number>
779 if (in_lagrange_product_form ==
true)
781 return lagrange_support_points.size();
786 return coefficients.size() - 1;
792 template <
typename number>
796 if (in_lagrange_product_form ==
false)
801 const unsigned int m = coefficients.size();
802 number value = coefficients.back();
803 for (
int k = m - 2; k >= 0; --k)
804 value = value * x + coefficients[k];
810 const unsigned int m = lagrange_support_points.size();
812 for (
unsigned int j = 0; j < m; ++j)
813 value *= x - lagrange_support_points[j];
814 value *= lagrange_weight;
821 template <
typename number>
822 template <
typename Number2>
825 const unsigned int n_derivatives,
829 if (in_lagrange_product_form ==
true)
834 const unsigned int n_supp = lagrange_support_points.size();
835 const number weight = lagrange_weight;
836 switch (n_derivatives)
840 for (
unsigned int d = 1;
d <= n_derivatives; ++
d)
842 for (
unsigned int i = 0; i < n_supp; ++i)
844 const Number2 v = x - lagrange_support_points[i];
852 for (
unsigned int k = n_derivatives; k > 0; --k)
853 values[k] = (values[k] * v + values[k - 1]);
865 number k_factorial = 1;
866 for (
unsigned int k = 0; k <= n_derivatives; ++k)
868 values[k] *= k_factorial * weight;
869 k_factorial *=
static_cast<number
>(k + 1);
882 for (
unsigned int i = 0; i < n_supp; ++i)
884 const Number2 v = x - lagrange_support_points[i];
887 values[0] = weight * value;
894 Number2 derivative = 0.;
895 for (
unsigned int i = 0; i < n_supp; ++i)
897 const Number2 v = x - lagrange_support_points[i];
898 derivative = derivative * v + value;
901 values[0] = weight * value;
902 values[1] = weight * derivative;
909 Number2 derivative = 0.;
911 for (
unsigned int i = 0; i < n_supp; ++i)
913 const Number2 v = x - lagrange_support_points[i];
914 second = second * v + derivative;
915 derivative = derivative * v + value;
918 values[0] = weight * value;
919 values[1] = weight * derivative;
920 values[2] =
static_cast<number
>(2) * weight * second;
931 const unsigned int m = coefficients.size();
932 std::vector<Number2> a(coefficients.size());
933 std::copy(coefficients.begin(), coefficients.end(), a.begin());
934 unsigned int j_factorial = 1;
939 const unsigned int min_valuessize_m =
std::min(n_derivatives + 1, m);
940 for (
unsigned int j = 0; j < min_valuessize_m; ++j)
942 for (
int k = m - 2; k >=
static_cast<int>(j); --k)
943 a[k] += x * a[k + 1];
944 values[j] =
static_cast<number
>(j_factorial) * a[j];
946 j_factorial *= j + 1;
950 for (
unsigned int j = min_valuessize_m; j <= n_derivatives; ++j)
956 template <
typename number>
957 template <
class Archive>
964 ar &in_lagrange_product_form;
965 ar &lagrange_support_points;
971 template <
typename Number>
978 Assert(alpha >= 0 && beta >= 0,
985 const Number xeval = Number(-1) + 2. * x;
991 p1 = ((alpha + beta + 2) * xeval + (alpha - beta)) / 2;
995 for (
unsigned int i = 1; i < degree; ++i)
997 const Number v = 2 * i + (alpha + beta);
998 const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
999 const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
1000 const Number a3 = v * (v + 1) * (v + 2);
1001 const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
1003 const Number pn = ((a2 + a3 * xeval) * p1 - a4 * p0) / a1;
1012 template <
typename Number>
1018 std::vector<Number> x(degree, 0.5);
1029 const Number tolerance =
1039 const unsigned int n_points = (alpha == beta ? degree / 2 : degree);
1040 for (
unsigned int k = 0; k < n_points; ++k)
1044 Number r = 0.5 - 0.5 * std::cos(static_cast<Number>(2 * k + 1) /
1047 r = (r + x[k - 1]) / 2;
1050 for (
unsigned int it = 1; it < 1000; ++it)
1053 for (
unsigned int i = 0; i < k; ++i)
1054 s += 1. / (r - x[i]);
1058 (alpha + beta + degree + 1) *
1063 const Number delta = f / (f * s - J_x);
1066 std::abs(delta) < tolerance)
1071 if (it == converged + 1)
1076 ExcMessage(
"Newton iteration for zero of Jacobi polynomial " 1077 "did not converge."));
1083 for (
unsigned int k = n_points; k < degree; ++k)
1084 x[k] = 1.0 - x[degree - k - 1];
static const unsigned int invalid_unsigned_int
void serialize(Archive &ar, const unsigned int version)
Polynomial< number > & operator*=(const double s)
void scale(const number factor)
void transform_into_standard_form()
bool operator==(const Polynomial< number > &p) const
Polynomial< number > derivative() const
void shift(const number2 offset)
unsigned int degree() const
Polynomial< number > & operator-=(const Polynomial< number > &p)
static std::vector< std::unique_ptr< const std::vector< double > > > recursive_coefficients
std::vector< Number > jacobi_polynomial_roots(const unsigned int degree, const int alpha, const int beta)
bool in_lagrange_product_form
number value(const number x) const
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
void print(std::ostream &out) const
Number jacobi_polynomial_value(const unsigned int degree, const int alpha, const int beta, const Number x)
static ::ExceptionBase & ExcMessage(std::string arg1)
static void multiply(std::vector< number > &coefficients, const number factor)
#define Assert(cond, exc)
Polynomial< number > & operator+=(const Polynomial< number > &p)
#define DEAL_II_NAMESPACE_CLOSE
virtual std::size_t memory_consumption() const
std::vector< number > coefficients
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Polynomial< number > primitive() const
std::vector< number > lagrange_support_points
static constexpr double PI
#define DEAL_II_NAMESPACE_OPEN
T min(const T &t, const MPI_Comm &mpi_communicator)
static ::ExceptionBase & ExcEmptyObject()
static ::ExceptionBase & ExcNotImplemented()
void copy(const T *begin, const T *end, U *dest)
T max(const T &t, const MPI_Comm &mpi_communicator)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)