37 template <
typename number>
40 , in_lagrange_product_form(
false)
46 template <
typename number>
48 : coefficients(n + 1, 0.)
49 , in_lagrange_product_form(
false)
55 template <
typename number>
58 : in_lagrange_product_form(
true)
65 for (
unsigned int i = 0; i <
supp.size(); ++i)
74 ExcMessage(
"Underflow in computation of Lagrange denominator."));
76 ExcMessage(
"Overflow in computation of Lagrange denominator."));
83 template <
typename number>
89 value(
x, values.size() - 1, values.data());
94 template <
typename number>
103 coefficients.resize(lagrange_support_points.size() + 1);
104 if (lagrange_support_points.empty())
105 coefficients[0] = 1.;
108 coefficients[0] = -lagrange_support_points[0];
109 coefficients[1] = 1.;
110 for (
unsigned int i = 1; i < lagrange_support_points.size(); ++i)
112 coefficients[i + 1] = 1.;
113 for (
unsigned int j = i;
j > 0; --
j)
114 coefficients[
j] = (-lagrange_support_points[i] * coefficients[
j] +
115 coefficients[
j - 1]);
116 coefficients[0] *= -lagrange_support_points[i];
119 for (
unsigned int i = 0; i < lagrange_support_points.size() + 1; ++i)
120 coefficients[i] *= lagrange_weight;
125 in_lagrange_product_form =
false;
126 lagrange_weight = 1.;
131 template <
typename number>
137 for (
typename std::vector<number>::iterator c = coefficients.begin();
138 c != coefficients.end();
148 template <
typename number>
155 if (in_lagrange_product_form ==
true)
157 number
inv_fact = number(1.) / factor;
159 for (
unsigned int i = 0; i < lagrange_support_points.size(); ++i)
161 lagrange_support_points[i] *=
inv_fact;
168 scale(coefficients, factor);
173 template <
typename number>
178 for (
typename std::vector<number>::iterator c = coefficients.begin();
179 c != coefficients.end();
186 template <
typename number>
190 if (in_lagrange_product_form ==
true)
191 lagrange_weight *= s;
194 for (
typename std::vector<number>::iterator c = coefficients.begin();
195 c != coefficients.end();
204 template <
typename number>
210 if (in_lagrange_product_form ==
true && p.in_lagrange_product_form ==
true)
212 lagrange_weight *= p.lagrange_weight;
213 lagrange_support_points.insert(lagrange_support_points.end(),
214 p.lagrange_support_points.begin(),
215 p.lagrange_support_points.end());
219 else if (in_lagrange_product_form ==
true)
220 transform_into_standard_form();
225 std::unique_ptr<Polynomial<number>>
q_data;
227 if (p.in_lagrange_product_form ==
true)
229 q_data = std::make_unique<Polynomial<number>>(p);
230 q_data->transform_into_standard_form();
237 unsigned int new_degree = this->degree() +
q->degree();
241 for (
unsigned int i = 0; i <
q->coefficients.size(); ++i)
242 for (
unsigned int j = 0;
j < this->coefficients.size(); ++
j)
251 template <
typename number>
263 if (in_lagrange_product_form ==
true)
264 transform_into_standard_form();
269 std::unique_ptr<Polynomial<number>>
q_data;
271 if (p.in_lagrange_product_form ==
true)
273 q_data = std::make_unique<Polynomial<number>>(p);
282 if (
q->coefficients.size() > coefficients.size())
283 coefficients.resize(
q->coefficients.size(), 0.);
285 for (
unsigned int i = 0; i <
q->coefficients.size(); ++i)
286 coefficients[i] +=
q->coefficients[i];
293 template <
typename number>
299 if (in_lagrange_product_form ==
true)
300 transform_into_standard_form();
305 std::unique_ptr<Polynomial<number>>
q_data;
307 if (p.in_lagrange_product_form ==
true)
309 q_data = std::make_unique<Polynomial<number>>(p);
310 q_data->transform_into_standard_form();
318 if (
q->coefficients.size() > coefficients.size())
319 coefficients.resize(
q->coefficients.size(), 0.);
321 for (
unsigned int i = 0; i <
q->coefficients.size(); ++i)
322 coefficients[i] -=
q->coefficients[i];
329 template <
typename number>
338 if (in_lagrange_product_form ==
true && p.in_lagrange_product_form ==
true)
339 return ((lagrange_weight == p.lagrange_weight) &&
340 (lagrange_support_points == p.lagrange_support_points));
341 else if (in_lagrange_product_form ==
true)
344 q.transform_into_standard_form();
345 return (
q.coefficients == p.coefficients);
347 else if (p.in_lagrange_product_form ==
true)
350 q.transform_into_standard_form();
351 return (
q.coefficients == coefficients);
354 return (p.coefficients == coefficients);
359 template <
typename number>
360 template <
typename number2>
379 const unsigned int n = d;
396 for (
unsigned int k = 0;
k < d; ++
k)
421 template <
typename number>
422 template <
typename number2>
429 if (in_lagrange_product_form ==
true)
431 for (
unsigned int i = 0; i < lagrange_support_points.size(); ++i)
432 lagrange_support_points[i] -= offset;
436 shift(coefficients, offset);
441 template <
typename number>
450 std::unique_ptr<Polynomial<number>>
q_data;
452 if (in_lagrange_product_form ==
true)
454 q_data = std::make_unique<Polynomial<number>>(*this);
455 q_data->transform_into_standard_form();
462 for (
unsigned int i = 1; i <
q->coefficients.size(); ++i)
470 template <
typename number>
476 std::unique_ptr<Polynomial<number>>
q_data;
478 if (in_lagrange_product_form ==
true)
480 q_data = std::make_unique<Polynomial<number>>(*this);
481 q_data->transform_into_standard_form();
489 for (
unsigned int i = 0; i <
q->coefficients.size(); ++i)
497 template <
typename number>
501 if (in_lagrange_product_form ==
true)
503 out << lagrange_weight;
504 for (
unsigned int i = 0; i < lagrange_support_points.size(); ++i)
505 out <<
" (x-" << lagrange_support_points[i] <<
")";
509 for (
int i = degree(); i >= 0; --i)
511 out << coefficients[i] <<
" x^" << i << std::endl;
516 template <
typename number>
530 template <
typename number>
534 std::vector<number>
result(n + 1, 0.);
541 template <
typename number>
548 template <
typename number>
549 std::vector<Polynomial<number>>
552 std::vector<Polynomial<number>> v;
553 v.reserve(degree + 1);
554 for (
unsigned int i = 0; i <= degree; ++i)
565 namespace LagrangeEquidistantImplementation
567 std::vector<Point<1>>
570 std::vector<Point<1>> points(n + 1);
572 for (
unsigned int k = 0;
k <= n; ++
k)
584 generate_equidistant_unit_points(n),
607 std::vector<double> &a)
611 unsigned int n_functions = n + 1;
613 const double *
x =
nullptr;
619 static const double x1[4] = {1.0, -1.0, 0.0, 1.0};
625 static const double x2[9] = {
626 1.0, -3.0, 2.0, 0.0, 4.0, -4.0, 0.0, -1.0, 2.0};
632 static const double x3[16] = {1.0,
656 for (
unsigned int i = 0; i < n_functions; ++i)
662 std::vector<Polynomial<double>>
667 return std::vector<Polynomial<double>>(
673 std::vector<Polynomial<double>> v;
674 for (
unsigned int i = 0; i <=
degree; ++i)
685 std::vector<Polynomial<double>>
688 std::vector<Polynomial<double>> p;
689 p.reserve(points.size());
691 for (
unsigned int i = 0; i < points.size(); ++i)
692 p.emplace_back(points, i);
714 for (
unsigned int i = 0; i <
k; ++i)
721 for (
unsigned int i = 0; i <
k; ++i)
728 std::vector<Polynomial<double>>
731 std::vector<Polynomial<double>> v;
733 for (
unsigned int i = 0; i <=
degree; ++i)
790 for (
unsigned int i = 2; i < p; ++i)
792 for (
unsigned int j = 0;
j < i - 1; ++
j)
795 for (
unsigned int j = 0;
j < i; ++
j)
806 for (
unsigned int j = 1;
j < i - 1; ++
j)
824 for (
int i = p; i > 0; --i)
833 std::vector<Polynomial<double>>
836 std::vector<Polynomial<double>> basis(p + 1);
838 for (
unsigned int i = 0; i <= p; ++i)
850 std::vector<std::unique_ptr<const std::vector<double>>>
924 std::vector<double>
c0(2);
928 std::vector<double>
c1(2);
935 std::make_unique<const std::vector<double>>(std::move(
c0));
937 std::make_unique<const std::vector<double>>(std::move(
c1));
945 std::vector<double>
c2(3);
954 std::make_unique<const std::vector<double>>(std::move(
c2));
970 std::vector<double>
ck(
k + 1);
976 for (
unsigned int i = 1; i <=
k - 1; ++i)
998 std::make_unique<const std::vector<double>>(std::move(
ck));
1004 const std::vector<double> &
1019 std::vector<Polynomial<double>>
1032 return std::vector<Polynomial<double>>(
1036 std::vector<Polynomial<double>> v;
1038 for (
unsigned int i = 0; i <=
degree; ++i)
1093 (*this) *= legendre;
1099 std::vector<Polynomial<double>>
1104 "degrees less than three"));
1105 std::vector<Polynomial<double>> basis(n + 1);
1107 for (
unsigned int i = 0; i <= n; ++i)
1144 for (
unsigned int q = 0;
q <
gauss.size(); ++
q)
1146 const double x =
gauss.point(
q)[0];
1148 for (
unsigned int j = 0;
j < degree - 3; ++
j)
1167 const unsigned int index)
1179 else if (degree == 1)
1192 else if (degree == 2)
1200 else if (index == 1)
1213 else if (degree == 3)
1230 else if (index == 1)
1240 else if (index == 2)
1247 else if (index == 3)
1291 for (
unsigned int m = 0; m <
degree - 3; ++m)
1299 else if (index == 1)
1302 for (
unsigned int m = 0; m <
degree - 3; ++m)
1315 std::vector<Point<1>> points(
degree);
1317 for (
unsigned int i = 0; i <
degree; ++i)
1335 else if (index >= 2 && index <
degree - 1)
1339 for (
unsigned int m = 0, c = 2; m <
degree - 3; ++m)
1349 else if (index ==
degree - 1)
1353 for (
unsigned int m = 0; m <
degree - 3; ++m)
1357 std::vector<Point<1>> points(
degree);
1359 for (
unsigned int i = 0; i <
degree; ++i)
1377 else if (index ==
degree)
1383 for (
unsigned int m = 0; m <
degree - 3; ++m)
1395 std::vector<Polynomial<double>>
1398 std::vector<Polynomial<double>> basis(
degree + 1);
1400 for (
unsigned int i = 0; i <=
degree; ++i)
HermiteInterpolation(const unsigned int p)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
HermiteLikeInterpolation(const unsigned int degree, const unsigned int index)
static std::vector< std::unique_ptr< const std::vector< double > > > recursive_coefficients
Hierarchical(const unsigned int p)
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
static void compute_coefficients(const unsigned int n, const unsigned int support_point, std::vector< double > &a)
LagrangeEquidistant(const unsigned int n, const unsigned int support_point)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Legendre(const unsigned int p)
std::vector< double > compute_coefficients(const unsigned int p)
Lobatto(const unsigned int p=0)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
static std::vector< Polynomial< number > > generate_complete_basis(const unsigned int degree)
Monomial(const unsigned int n, const double coefficient=1.)
static std::vector< number > make_vector(unsigned int n, const double coefficient)
number value(const number x) const
bool operator==(const Polynomial< number > &p) const
std::vector< double > coefficients
Polynomial< number > primitive() const
Polynomial< number > & operator+=(const Polynomial< number > &p)
Polynomial< number > derivative() const
void transform_into_standard_form()
void scale(const number factor)
Polynomial< number > & operator-=(const Polynomial< number > &p)
std::vector< number > lagrange_support_points
void shift(const number2 offset)
void print(std::ostream &out) const
bool in_lagrange_product_form
static void multiply(std::vector< number > &coefficients, const number factor)
Polynomial< number > & operator*=(const double s)
virtual std::size_t memory_consumption() const
unsigned int degree() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcZero()
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcEmptyObject()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_ASSERT_UNREACHABLE()
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
std::vector< Point< 1 > > generate_equidistant_unit_points(const unsigned int n)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 > > &points)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)